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RECENT  ADVANCES  IN  ACHIEVING  TEXTBOOK  MULTIGRID  EFFICIENCY  FOR 
COMPUTATIONAL  FLUID  DYNAMICS  SIMULATIONS 

ACHI  BRANDT*,  BORIS  DISKINt,  AND  JAMES  L.  THOMAS^ 

Abstract.  Recent  advances  in  achieving  textbook  multigrid  efficiency  for  fluid  simulations  are  presented. 
Textbook  multigrid  efficiency  is  defined  as  attaining  the  solution  to  the  governing  system  of  equations  in 
a  computational  work  which  is  a  small  multiple  of  the  operation  counts  associated  with  discretizing  the 
system.  Strategies  are  reviewed  to  attain  this  efficiency  by  exploiting  the  factorizability  properties  inherent 
to  a  range  of  fluid  simulations,  including  the  compressible  Navier-Stokes  equations.  Factorizability  is  used  to 
separate  the  elliptic  and  hyperbolic  factors  contributing  to  the  target  system;  each  of  the  factors  can  then  be 
treated  individually  and  optimally.  Boundary  regions  and  discontinuities  are  addressed  with  separate  (local) 
treatments.  New  formulations  and  recent  calculations  demonstrating  the  attainment  of  textbook  efficiency 
for  aerodynamic  simulations  are  shown. 

Key  words,  textbook  multigrid  efficiency,  factorizable  systems  of  differential  equations,  principal  lin¬ 
earization,  factorizable  discretizations,  distributed  relaxation 

Subject  classification.  Applied  and  Numerical  Mathematics 

1.  Introduction.  Considerable  progress  over  the  past  thirty  years  has  been  made  in  the  development 
of  large-scale  computational  fluid  dynamics  (CFD)  solvers  for  the  Euler  and  Navier-Stokes  equations.  Com¬ 
putations  are  used  routinely  to  design  the  cruise  shapes  of  transport  aircraft  through  complex-geometry 
simulations  involving  the  solution  of  25-100  million  equations;  in  this  arena,  the  number  of  wind-tunnel  tests 
for  a  new  design  has  been  substantially  reduced.  However,  simulations  of  the  entire  flight  envelope  of  the 
vehicle,  including  maximum  lift,  buffet  onset,  flutter,  and  control  effectiveness,  have  not  been  as  successful 
in  eliminating  the  reliance  on  wind-tunnel  testing.  These  simulations  involve  unsteady  flows  with  more 
separation  and  stronger  shock  waves  than  at  cruise.  The  main  reasons  limiting  further  inroads  of  CFD  into 
the  design  process  are:  (1)  the  reliability  of  turbulence  models  and  (2)  the  time  and  expense  of  the  numer¬ 
ical  simulation.  Because  of  the  prohibitive  resolution  requirements  of  direct  simulations  at  high  Reynolds 
numbers,  transition  and  turbulence  modeling  is  expected  to  remain  an  issue  for  the  near  term  [41].  The 
focus  of  this  paper  addresses  the  latter  problem  by  attempting  to  attain  optimal  efficiencies  in  solving  the 
governing  equations.  Typically  current  CFD  codes  based  on  the  use  of  multigrid  acceleration  techniques  and 
multistage  Runge-Kutta  time-stepping  schemes  are  able  to  converge  lift  and  drag  values  for  cruise  configura¬ 
tions  within  approximately  1000  residual  evaluations.  More  complexity  in  the  geometry  or  physics  generally 
requires  many  more  residual  evaluations  to  converge,  and  sometimes  convergence  cannot  be  attained.  An 
optimally  convergent  method  is  defined  [5,  6,  8,  9]  as  having  textbook  multigrid  efficiency  (TME),  meaning 
the  solutions  to  the  governing  system  of  equations  are  attained  in  a  computational  work  which  is  a  small 
(less  than  10)  multiple  of  the  operation  count  in  the  discretized  system  of  equations  (residual  evaluations). 
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Thus,  there  is  a  potential  gain  of  more  than  two  orders  of  magnitude  in  operation  count  reduction  if  TME 
could  be  attained.  For  staggered-grid  formulations  of  incompressible  flow  equations,  robust  and  relatively 
efficient  multigrid  solvers  have  already  been  developed  [30,  31,  32];  their  efficiency  is  still  about  an  order  of 
magnitude  behind  TME. 

Because  the  governing  equations  are  a  set  of  coupled  nonlinear  conservation  equations  with  discontinuities 
(shocks,  slip  lines,  etc.)  and  singularities  (flow-  or  grid-induced),  the  difficulties  are  numerous.  The  TME 
methodology  insists  that  each  of  the  difficulties  should  be  isolated,  analyzed,  and  solved  systematically 
using  a  carefully  constructed  series  of  model  problems.  To  be  efficient,  multigrid  solvers  for  general  systems 
of  partial  differential  equations  must  adequately  address  three  types  of  errors:  (1)  high-frequency  error 
components,  (2)  uniformly  smooth  error  components,  and  (3)  characteristic  error  components.  The  latter 
are  (usually  smooth)  error  components  that  are  much  smoother  in  characteristic  directions  than  in  other 
directions.  Standard  multigrid  methods  that  are  efficient  for  elliptic  problems  recognize  and  separate  the 
treatment  of  high-frequency  and  smooth  error  components.  The  former  are  efficiently  reduced  in  relaxation; 
the  latter  are  well  approximated  on  coarse  grids  and,  hence,  eliminated  through  the  coarse-grid  correction. 
The  efficiency  of  classical  multigrid  methods  severely  degrades  for  nonelliptic  problems  because  characteristic 
components  cannot  be  adequately  approximated  on  coarse  grids  [4,  11,  15]. 

If  the  target  discretization  is  /i-elliptic  (or  semi-/i-elliptic) ,  the  high  frequency  error  components  can 
still  be  reduced  by  a  local  (or  block- wise)  relaxation  procedure.  By  definition  [4,  55],  a  discrete  scalar 
(not  necessarily  elliptic)  operator  L[u]  possesses  a  good  measure  of  /i-ellipticity  if  the  absolute  value  of  its 
symbol  |T(^)|  =  is  well  separated  from  zero  for  all  high-frequency  Fourier  modes.  Here 

j  =  Ux,jy,jz)  are  the  grid  indexes  and  0  =  {9x,0y,9z),0  <  J^a,],  J^j,],  J^z]  <  tt  are  normalized  Fourier 
frequencies.  High-frequency  Fourier  modes  are  the  modes  satisfying  max(|0a;j,  \9y\,  j^zj)  >  f  •  For  systems, 
the  measure  of  /i-ellipticity  is  defined  as  the  absolute  value  of  the  determinant  of  the  operator  matrix. 

Standard  coarse-grid  corrections  are  efficient  for  uniformly  smooth  error  components,  even  for  nonelliptic 
problems.  An  effective  reduction  of  characteristic  error  components  can  be  achieved  either  by  designing  a 
proper  relaxation  scheme  reducing  not  only  high-frequency  but  smooth  error  components  as  well  (e.g.,  by 
downstream  ordering  of  relaxation  steps)  or  by  adjusting  coarse-grid  operators  for  a  better  characteristic- 
component  approximation. 

Multigrid  methods  efficiently  reducing  all  the  three  aforementioned  types  of  error  have  been  developed  for 
scalar  nonelliptic  operators  [11,  12,  20,  21,  22].  TME  for  systems  of  equations  can  be  attained  by  exploiting 
the  factorizability  property  of  the  governing  equations.  The  factorizability  of  the  Navier-Stokes  equations  is 
manifested  by  the  fact  that  the  determinant  of  the  matrix  of  the  differential  operators  consists  of  separable 
factors.  Exploiting  the  factorizability  property  in  discrete  computations  reduces  the  problem  of  relaxing  a 
complicated  system  of  discretized  coupled  differential  equations  to  relaxation  of  simpler  factors  constituting 
the  system  determinant.  This  approach  is  quite  distinct  from  most  approaches  to  accelerate  convergence 
because  for  steady-state  flows,  the  factors  are  treated  directly  rather  than  through  pseudo-time  marching 
methods.  Time-dependent  flow  solvers  can  be  constructed  within  this  approach  as  well  and  in  principle 
are  simpler  to  develop  than  steady-state  solvers.  A  list  of  envisioned  difficulties  and  possible  solutions  in 
attaining  TME  for  CFD  simulations  is  discussed  elsewhere  [9]. 

This  paper  is  organized  as  follows.  The  foundations  of  the  methodology  for  attaining  TME  in  CFD 
simulations  are  discussed  in  Section  2,  including  the  concept  of  principal  linearization  and  illustrations  of 
the  factorizability  of  various  fluid  dynamic  equations.  Two  strategies  for  exploiting  the  factorizability  are 
presented  in  the  next  two  sections.  Reformulation  of  the  differential  equations  is  discussed  in  Section  3.  An 
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alternative,  more  general,  distributed  relaxation  approach  is  discussed  in  Section  4.  Efficient  methods  for 
relaxing  scalar  factors  are  presented  in  Section  5.  Section  6  compares  and  evaluates  the  available  analytical 
tools.  Section  7  summarizes  recent  advances  in  formulation  and  demonstration  TME.  Concluding  remarks 
are  given  in  the  final  Section  8. 

2.  Foundations  for  Textbook  Multigrid  Efficiency.  The  basic  framework  for  TME  solvers  is  full 
multigrid  (EMG)  algorithms  [5,  6,  17,  46,  55,  56].  In  EMG  algorithms,  the  solution  process  is  started  on 
a  very  coarse  grid  where  the  computational  cost  of  solution  is  negligible.  The  coarse-grid  solution  is  then 
interpolated  to  the  next  fine  grid  to  form  an  initial  approximation.  Eew  multigrid  full  approximation  scheme 
(EAS)  cycles,  or  possibly  just  one,  are  performed  next  to  obtain  an  improved  fine-grid  solution  approximation. 
Then,  the  process  proceeds  to  finer  grids  until  the  solution  on  the  target  finest  grid  is  achieved. 

In  the  solution  of  highly  nonlinear  problems,  a  good  initial  guess  is  important.  A  general  way  to  obtain 
such  an  initial  guess  is  by  continuation,  in  which  the  solution  to  the  target  problem  is  approached  through 
the  solutions  of  a  sequence  of  parameterized  problems.  Usually  the  problem  starting  the  continuation  process 
is  easy  to  solve,  and  the  difficulty  gradually  increases  as  the  control  parameter  approaches  the  target  value; 
this  continuation  process  can  often  be  integrated  into  an  EMG  solver.  Eor  example,  with  viscosity  as  the 
control  parameter,  at  the  coarse  grids  more  artificial  viscosity  can  be  used,  then  gradually  be  taken  out  as 
the  algorithm  proceeds  to  finer  levels.  Such  EMG  continuation  is  often  natural  because  larger  numerical 
viscosity  is  introduced  on  coarse  grids,  even  without  aiming  at  continuation. 

A  version  named  A-EMG  algorithm  provides  the  device  needed  for  optimal  adaptive  local  refinement. 
Efficient  multigrid  solvers  based  on  this  approach  have  been  demonstrated  [2] . 

The  objective  of  EMG  algorithms  (and  TME  methods  in  particular)  is  fast  convergence  to  the  solution 
of  the  differential  equations,  not  necessarily  fast  asymptotic  residual  convergence.  The  natural  solution 
tolerance  is  the  discretization  error  defined  as  the  difference  between  the  exact  solutions  of  discrete  and 
differential  problems.  Thus,  the  quality  of  a  solution  approximation  on  a  given  grid  can  be  measured  by  the 
relative  magnitude  of  algebraic  errors  in  comparison  with  the  discretization  error  level.  The  algebraic  error 
is  defined  as  the  difference  between  the  exact  and  approximate  solutions  of  the  discrete  problem.  On  any 
grid  in  an  EMG  algorithm,  we  expect  the  algebraic  errors  after  few  multigrid  cycles  to  be  always  less  than 
the  discretization  error. 

On  the  other  hand,  a  fast  residual  convergence  is  considered  as  an  important  monitoring  tool.  In  many 
practical  cases,  it  is  possible  to  develop  a  solver  exhibiting  fast  residual  convergence  rates  without  compro¬ 
mising  TME.  Note  however  that  sometimes  the  quality  of  the  target-grid  solution  can  be  much  improved 
by  double  discretization  methods.  A  double  discretization  method  applies  for  relaxation  a  discretization 
scheme  that  is  different  from  the  scheme  used  for  calculating  residuals  transferred  to  the  coarse  grid.  Zero 
target-grid  residuals  might  not  be  the  aim  in  this  case. 

The  goal  of  this  paper  is  to  review  solution  strategies  leading  to  TME  solution  of  fluid  mechanics 
equations.  The  most  general  system  we  consider  here  is  the  time-dependent  compressible  Navier-Stokes 
equations  written  as 

(2.1)  aiQ  +  R(Q)  =  0, 

where  the  conserved  variables  are  Q  =  {pu,  pv,  pw,  p,  pE)^ ,  representing  the  momentum  vector,  density,  and 
total  energy  per  unit  volume,  and 

(2.2)  R(Q)  =  a,F(Q)+a,G(Q)+a,H(Q), 
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where 


F(Q)  = 


G(Q)  = 


H(Q)  = 


fm?  +  p  —  2iidxU  —  A(v  •  u) 

puv  —  p(dxV  +  dyu) 
puw  —  p(dxW  +  dzu) 
pu 

y  puE  +  up  —  Au(v  •  u)  —  pTi  —  ndxf. 

/  puv  —  p{dxV  +  dyu) 

pv"^  +  p  -  2pdyV  —  A(v  •  u) 

pvw  —  p(dyW  +  dzv) 
pv 

\  pvE  +  vp  —  Xv(sj  •  u)  —  pT2  —  Kdye 

/  puw  —  p{dxW  +  dzu) 

pvw  —  p(dyW  +  dzv) 

pw"^  +  p  —  2pdzW  —  A( V  •  u) 

pw 

\  pwE  +  wp  —  Aw(v  •  u)  —  pTs  —  Kdze 


\ 

\ 

\ 

y 


Ti  =  2udxU  +  v(dxV  +  dyu)  +  w(dxW  +  dzu), 
T2  =  2vdyV  +  u(dxV  +  dyu)  +  w(dyW  +  dzv), 
Ts  =  2wdzW  +  u(dxW  +  dzu)  +  v(dyW  +  dzv), 


p  and  A  are  viscosity  coefficients,  and  k  is  the  coefficient  of  heat  conductivity. 

A  basic  step  in  developing  an  efficient  multigrid  algorithm  is  to  design  an  efficient  relaxation  procedure. 
For  nonlinear  problems,  the  relaxation  updates  to  a  current  solution  approximation  are  usually  computed 
through  Newton  iterations.  The  full  Newton  linearization  of  the  Navier-Stokes  equations  (2.1)  is  a  very 
complicated  operator,  and  its  solution  (inversion)  is  too  costly  for  practical  applications.  To  reduce  the  com¬ 
putational  cost  without  compromising  efficiency,  one  can  opt  for  relaxation  of  a  principal  linearization.  The 
principal  linearization  of  a  scalar  equation  contains  the  linearization  terms  that  make  a  major  contribution 
to  the  residual  per  unit  change  in  the  unknown  variable.  The  principal  terms  thus  generally  depend  on  the 
scale,  or  mesh  size,  of  interest.  For  example,  the  discretized  highest  derivative  terms  are  principal  on  grids 
with  small  enough  mesh  size.  For  a  discretized  system  of  differential  equations,  the  principal  terms  are  terms 
that  contribute  to  the  principal  terms  of  the  system  determinant. 

To  illustrate  the  idea  of  principal  linearization,  consider  a  nonlinear  discrete  thin-layer  approximation 
of  the  convection-diffusion  operator,  in  which  the  flow  is  parallel  to  the  boundary  (a;-direction)  and  only  the 
viscous  terms  associated  with  variations  in  the  y-coordinate  normal  to  the  boundary  are  retained 


(2.3)  N{u)  =  udxU  -  vdyyU, 

where  3^  and  dyy  are  discrete  approximations  to  the  first  a;-directional  derivative  and  to  the  second  y- 
directional  derivative,  respectively.  A  full  Newton  linearization  (assuming  constant  viscosity)  for  a  correction 
Su  has  three  terms 

(2.4)  ^  udxSu  +  {dxu)Su  —  vdyySu. 
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To  evaluate  principality  of  the  terms,  one  can  start  from  an  exact  discrete  solution;  for  example,  the  function 
Su  =  0  is  the  exact  solution  of  the  homogeneous  equation  ^Su  =  0.  A  unit  change  in  the  unknown  variable, 
Su,  is  defined  as  a  perturbation  of  the  solution  value  at  one  grid  point.  Introducing  this  perturbation,  the 
residual  function  becomes  nonzero  in  the  vicinity  of  the  perturbed  point.  One  can  directly  check  which  of  the 
terms  of  the  full  linearization  operator  make  major  residual  contributions.  Three  situations  are  encountered: 

uh^ 

•  High  cell  Reynolds  number  {jrif'  ^  !)•  If  velocity  function  u  is  smooth  and  non-degenerate, 
i.e.,  the  magnitude  of  velocity  deviations  in  neighboring  grid  points  is  less  than  the  local  velocity 
magnitude,  then  the  major  contribution  to  the  residual  function  is  0{u/hx)  and  comes  from  the 
term  ud^Su]  the  second  term  contribution  is  0(3^u);  and  the  viscous  term  contributes  0{i'/hy) 
which  is  much  less  than  0{u/hx)-  Thus,  only  the  first  term  ud^Su  is  principal.  However,  if  the 
velocity  field  is  not  smooth,  either  because  of  a  coarse  mesh  or  proximity  to  a  discontinuity,  or  if  the 
absolute  velocity  value  is  small  (stagnation  flows),  then  the  second  term  becomes  principal  as  well. 


•  Low  cell  Reynolds  number  1).  The  major  residual  contribution  comes  from  the  viscous  term 

which  is  the  only  principal  term. 

uh^ 

•  Medium  cell  Reynolds  number  =  0(1))-  This  situation  corresponds  to  the  usual  boundary 
layer  assumption  when  convection  balances  diffusion.  For  smooth  non-degenerate  flows,  the  only 
subprincipal  term  is  the  second;  the  first  and  the  third  terms  are  principal.  For  nonsmooth  or 
degenerate  flows,  all  three  terms  are  principal. 

As  an  example  of  a  system  of  nonlinear  flow  equations,  we  consider  a  discrete  operator  corresponding  to 
a  one-dimensional  steady-state  compressible  inviscid  flow: 


(2.5) 


iV(q) 


I'  ud^u  + 

'ypdxU  +  udxP 
V  (7  -  l)e3(lu  +  ud^e 


where  q  =  (u,p,  e)^  represents  velocity,  pressure,  and  internal  energy  and  7  is  the  ratio  of  specific  heats. 
The  full  Newton  linearization  of  this  operator  is  given  by 


(2.6) 


(d^u)  + ud^ 
ipdx  + 

(7  -  l)ea^  -b  (a^e) 


(7  -  i)<ia^  -  ifi) 

7(a(fu)  -b  ud^ 

0 


(7-1)^ 

0 

(7-  i)(3^u)  -bua^ 


f  Su 
Sp 

\  Se  ) 


Assuming  a  smooth  non-degenerate  solution  q,  the  first  simplification  step  is  to  eliminate  lower  derivative 
terms  from  each  entry  of  the  matrix  (2.6).  This  simplification  leads  to  an  approximate  linearization  as 


(2.7) 


f  ud^ 

(7-i)fa,^ 

(7-l)^t^ 

\ 

/  Su  \ 

ipdx 

udx 

0 

1  1 

\  (7  -  i)ea^ 

0 

udx 

) 

\  Se  j 

The  determinant  of  the  matrix  (2.7)  is 

^/l 


udiiu^didi  -  -  (7  -  lY-{d'^xP)d'i), 


where  the  sound  speed  c  relates  to  the  internal  energy  e  as  =  7(7  —  l)e.  Because  p  is  non-degenerate,  the 
last  term  in  the  parentheses  is  subprincipal  in  comparison  to  the  other  two  terms.  The  third  element  in  the 
first  row  of  the  matrix  (2.7)  does  not  contribute  to  the  principal  part  of  the  determinant  operator;  therefore 
the  principal  linearization  is  defined  as 
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(2.8) 


/  ud^ 

(7-i)fa,^ 

0 

\ 

^  5u  ^ 

ipdt 

ud^ 

0 

Sp 

V  (7  -  i)ea^ 

0 

ud^ 

/ 

V  J 

The  notion  of  principal  linearization  is  essentially  based  on  the  discrete  formulation.  The  principal  part 
of  a  differential  operator  may  be  defined  as  the  limit  of  the  principal  part  of  the  corresponding  discrete 
operator  as  the  mesh  size  h  tends  to  zero.  So  for  smooth  non-degenerate  flows,  the  principal  terms  of  the 
differential  equations  are  the  highest  derivatives. 

TME  for  solution  of  the  Navier-Stokes  system  of  differential  equations  can  be  achieved  by  exploiting  the 
system  factorizability.  To  illustrate  the  factorizability  property,  examples  are  given  below  for  various  fluid 
mechanics  regimes.  In  all  cases,  we  assume  a  smooth  non-degenerate  solution  as  defined  previously. 

Incompressible  Navier-Stokes  Equations.  The  steady-state  incompressible  Navier-Stokes  equations  can 
be  written  as 


(2.9) 


Qi,u  -I-  VP  =  0, 

V  u  =  0, 


where  u  =  {u,v,w)'^  is  the  velocity  vector  and  Q„  =  u  ■  —  vA  is  a  convection-diffusion  operator.  The 

principal  linearization  operator  is  given  by 


(2.10) 


where  velocity  u  is  fixed  in  the  linearized  convection-diffusion  operator  Q^. 
(2.11)  AetE  =  -QlA. 


^  Su 

0  0  dx 

^  Su 

Sv 

0  0  dy 

Sv 

Sw 

0  0  Q,  d. 

Sw 

V  ) 

_  dx  dy  0  _ 

[  Sp  ) 

Compressible  Euler  Equations.  A  nonconservative  form  of  the  Euler  equations  is  given  by 

Qu-l-i  VP  =  0, 
pc^  V  'U  +  Qp  =  0, 

Y  V  'U  +  Qe  =  0, 

where  Q  =  Qo  denotes  the  particular  case  of  Q„  with  zero  {v  =  0)  physical  diffusion,  and  density,  p,  pressure, 
p,  sound  speed,  c,  and  internal  energy,  e,  are  related  as 


(2.12) 

(2.13) 


P  =  (7  -  l)pe, 

=  jp/p. 


The  principal  linearization  is  given  by 


(2.14) 


/ 

Su 

] 

Q 

0 

0 

lid. 

0 

f 

Su 

\ 

Sv 

1 

0 

Q 

0 

idy 

0 

Sv 

Sw 

= 

0 

0 

Q 

id^ 

0 

Sw 

Sp 

1 

pc^dx 

pc^dy 

pc^d^ 

Q 

0 

Sp 

V 

Se 

J 

—d 

L  'y 

—d 

7 

—  d 

0 

Q  _ 

V 

Se 

/ 
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The  determinant  of  the  matrix  operator  L  is 


(2.15)  detL  =  [q2  ^ 

where  A  is  the  Laplace  operator,  and  —  c^A  represents  the  full-potential  operator. 

Compressible  Navier-Stokes  Equations.  The  nonconservative  formulation  corresponding  to  the  steady- 
state  version  of  (2.1)  is  given  by 

((u  •  v)  -  f  A  -  -  ^{d^yv  +  d^^w)  +  =  0, 

((u  •  V)  -  f  A  -  ^dyy^v  -  ^{d^yu  +  dy^w)  +  jdyp  =  0, 

((u  •  V)  -  f  A  -  -  ^(d^zu  +  dyzv)  +  jd^p  =  0, 

pc^iV  •  u)  -I-  (u  •  S7)p  +  (7  -  l)(-/tAe  -I-  $)  =  0, 

^(V  •  u)  -I-  (u  •  v)p  -  f  Ae  -I-  =  0, 

where 

$  =  p(2{dxuy  +  2{dyvy  +  2{dzwy  +  {dxV  +  dyu)"^  +  {dxW  +  dzu)"^  +  [dyW  +  dzv)^'j  +  X{dxU  +  dyV  +  dzwy . 

Assuming  constant  viscosity  and  heat  conduction  coefficients,  the  principal  linearization  operator  L, 
keeping  the  terms  principal  on  both  the  viscous  and  inviscid  scales,  is  given  by 


QiL  —  ^dxx 

p  ^ 

p  ^xy 

--9 

pUxz 

\dx 

0 

--d 

p  ^xy 

—  kd 

p^y^ 

-pdy 

0 

L  = 

—  kf) 

p^xz 

—  kf) 

pOyz 

QiL  -  Xdzz 

P  P 

0 

pc^dx 

pc  dy 

pc^dz 

Q 

(1  —  7)kA 

—d 

L  7 

—d 

-dz 

7 

0 

^7  - 

(2.17) 


det  L  =  Ql  \—A^  +  Qi-c^A  +  -  Q 

p  L  7p 


2^  +  X  +  p 


A  +  Q' 


where,  nondimensionalizing  by  density  and  sound  speed  and  applying  Stokes  hypothesis  for  the  bulk  viscosity 
term,  the  coefficients  become  pj p  =  M^o/ip  Re),/t  =  Moo7/(Re  Pr),  and  X  =  X  +  p  =  p/S;  M^o  is  the  free 
stream  Mach  number,  and  Re  and  Pr  are  Reynolds  and  Prandtl  numbers  respectively. 

The  approaches  exploiting  the  factorizability  property  for  efficient  solution  of  the  Navier-Stokes  equa¬ 
tions  may  be  divided  into  two  categories:  (1)  reformulating  the  target  differential  equations  so  that  the 
principal  linearization  of  the  new  formulation  becomes  uncoupled  (usually  triangular  with  the  factors  of  the 
determinant  on  the  main  diagonal)  and  (2)  modifying  the  equations  for  computing  relaxation  updates  while 
keeping  the  original  formulation  for  computing  residuals. 

For  the  subsonic  compressible  Euler  equations,  the  first  TME  solvers  exploiting  factorizability  of  the 
system  have  been  developed  by  Ta’asan  [48,  49,  50].  These  solvers  represent  examples  of  the  reformulation 
approach.  New  canonical  variables  have  been  introduced,  and  in  these  variables,  the  Euler  system  of  equa¬ 
tions  becomes  block  upper  triangular  with  the  main  diagonal  blocks  consisting  of  the  basic  components  of 
the  system.  Another  reformulation  approach  toward  achieving  TME  for  solution  of  the  Euler  and  incom¬ 
pressible  Navier-Stokes  equations  is  based  on  the  pressure-equation  formulation  which  effectively  separates 
elliptic  and  hyperbolic  factors  of  the  system  [36,  43,  45]. 
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The  approaches  from  the  second  category  are  more  general,  allowing  considerable  freedom  in  relaxation 
scheme  design  because  different  schemes  may  be  applied  to  different  flow  regions.  However,  the  design  itself 
is  relatively  simple  only  if  the  target  discretization  scheme  is  also  factorizable,  i.e.,  the  determinant  of  the 
discrete  principal  linearization  can  be  represented  as  a  product  of  discrete  factors,  each  of  them  approximating 
a  corresponding  factor  of  the  determinant  of  the  differential  equations.  A  stumbling  block  that  has  prevented 
fast  progress  in  developing  TME  solvers  was  the  lack  of  factorizable  discretizations  for  many  important 
application  areas  in  fluid  mechanics.  Among  the  widely  known  discretization  schemes,  only  staggered- 
grid  formulations  for  incompressible  and  subsonic  compressible  flow  regimes  are  conveniently  factorizable. 
Some  centrally-differenced  collocated-grid  formulations  are  factorizable  as  well,  but  the  factors  obtained 
in  the  corresponding  discrete  determinant  are  not  often  easily  treated.  The  search  for  new  factorizable 
discretization  schemes  (see  Summary  of  Recent  Progress  below)  is  chiefly  motivated  by  the  need  to  derive 
discrete  schemes  with  the  resulting  discretizations  of  scalar  factors  satisfying  some  desired  properties  (e.g., 
stability,  correct  alignment  with  the  physical  anisotropies,  compactness,  availability  of  an  efficient  relaxation 
scheme,  etc.)  Development  of  suitable  factorizable  discrete  schemes  for  the  Navier-Stokes  equation  is  a 
challenging  task.  Much  of  the  recent  progress  in  achieving  TME  for  CED  simulations  is  because  new  families 
of  general  factorizable  collocated-grid  discretization  schemes  are  emerging  [26,  35,  38,  42,  44].  The  next  two 
sections  present  some  details  of  methods  from  the  two  categories. 


3.  Equation  Reformulation  Strategies.  As  mentioned  previously,  the  first  TME  solvers  exploiting 
factorizability  of  the  system  have  been  developed  in  [48,  49,  50].  The  original  equations  were  reformulated 
in  terms  of  canonical  forms,  in  which  the  subsystems  governed  by  hyperbolic  operators  are  distinguished 
and  treated  separately,  both  in  discretization  and  relaxation,  from  those  governed  by  elliptic  operators.  The 
canonical  variables  for  two  dimensions  are  velocity  {u,v),  entropy  s,  and  total  enthalpy  H.  The  elliptic 
operators  are  discretized  with  /i— elliptic  centered  differences  and  solved  with  point  relaxation  and  coarse 
grid  corrections;  the  hyperbolic  operators  are  discretized  with  upwind  schemes  and  solved  by  marching 
techniques.  Ta’asan  [49]  was  able  to  demonstrate  solutions  for  the  subsonic  compressible  Euler  equations 
which  converged  with  the  same  rates  as  the  solution  of  the  scalar  full  potential  equation.  An  additional 
advantage  shown  for  this  formulation  was  that  the  total  artificial  viscosity  error  was  smaller  than  with  other 
schemes  because  the  upwinding  was  only  used  for  the  hyperbolic  subsystems.  The  main  disadvantage  of 
this  formulation  is  that  it  is  not  easily  generalized  to  viscous  and  unsteady  problems,  especially  in  three 
dimensions. 

An  alternative  pressure-equation  formulation  [45]  for  the  incompressible  Navier-Stokes  equations  (2.9) 
effectively  separates  the  elliptic  and  hyperbolic  factors  of  the  system.  The  continuity  equation  is  replaced 
with  an  equation  for  the  pressure,  as 


(3.1) 


v((3i^u  -I-  S7p)  -  Qi^iv  •  u)  =  0. 


Assuming  a  smooth  non-degenerate  flow,  the  principal  linearization  taken  in  the  limit  as  mesh  size  h  tends 
to  zero  is  an  upper  triangular  matrix  with  the  main  diagonal  composed  of  the  linearized  convection-diffusion 
and  Laplace  operators. 


(3.2) 


Qi,  0  0  dx 

0  0  dy 

0  0  Q,  d, 

0  0  0  A 

The  relaxation  scheme  is  defined  as 


(3.3)  L(5q=-i?(q), 


where  i?(q)  is  the  new  nonlinear  formulation  of  the  incompressible  Navier-Stokes  equations. 

The  determinant  of  the  reformulated  system  is  QlA  as  compared  to  QlA  of  the  original  system.  Thus 
additional  boundary  conditions  which  enforce  zero-divergence  need  to  be  applied.  The  equations  are  uncou¬ 
pled  everywhere  except  at  the  boundaries  and  some  local  relaxation  is  needed  to  relax  the  equations  in  this 
region.  Some  two-dimensional  results  are  shown  subsequently  demonstrating  the  efficiency  of  this  approach. 
Although  unexploited  as  yet,  the  approach  applies  equally  well  to  time-dependent  flows.  This  approach  has 
met  difficulties  in  generalizing  to  viscous  compressible  flows. 

4.  Distributed  Relaxation  Strategies.  The  most  general  procedure  exploiting  factorizability  of  the 
target  Navier-Stokes  equations  is  distributed  relaxation.  The  general  framework,  first  introduced  in  [54]  can 
be  outlined  as  follows. 

In  general,  the  simplest  form  of  the  differential  Navier-Stokes  equations  corresponds  to  nonconservative 
equations  expressed  in  primitive  variables,  e.g.,  taken  as  the  set  composed  of  velocity,  pressure,  and  internal 
energy,  q  =  (u,v,w,p,  e)^.  For  a  perfect  gas,  the  primitive  and  conservative  variables  are  connected  through 
(2.12),  (2.13),  and 


(4.1)  e  =  E  -  +v^ +w^'j. 

The  time-dependent  nonconservative  equations  are  found  readily  by  transforming  the  time-dependent  con¬ 
servative  equations. 

|§[aiQ  +  R]  =  0, 
a^q  +  |§R  =0, 

where  is  the  Jacobian  matrix  of  the  transformation.  For  steady-state  equations,  the  time  derivative  is 
dropped.  In  an  iterative  procedure,  the  correction  (5q  =  q"+^  —  q",  where  n  is  an  iteration  counter,  can  be 
computed  from  the  equation 


(4.2)  = 

where  the  right  side  of  (4.2)  is  a  linear  combination  of  the  conservative  residuals,  and  L  is  the  principal 
linearization  of  the  nonconservative  operator  at  the  scale  h. 

While  significantly  simplified  by  retaining  only  principal  terms,  the  system  (4.2)  is  still  a  set  of  coupled 
equations  containing  elliptic  and  hyperbolic  components.  Therefore,  collective  Gauss-Seidel  relaxation  of  L 
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is  not  often  effective.  The  distributed  relaxation  method  replaces  (5q  in  (4.2)  by  M(5w. 

(4.3)  LM(5w  = 

The  resulting  matrix  L  M  becomes  lower  triangular.  The  diagonal  elements  of  L  M  are  composed  ideally  of 
the  separable  factors  of  the  matrix  L  determinant.  These  factors  are  scalar  differential  operators  of  first  or 
second  order,  so  their  efficient  relaxation  is  a  much  simpler  task  than  relaxing  the  entire  system  associated 
with  L.  In  relaxing  scalar  factors,  the  changes  introduced  in  the  “ghost”  variables  Sw  (the  variables  (5w 
are  “ghost”  because  they  need  not  be  explicitly  used  in  computations)  during  relaxation  are  distributed, 
with  the  pattern  of  the  distribution  matrix  M,  to  the  primitive  variables.  To  obtain  the  optimal  (textbook) 
efficiency,  relaxation  of  each  factor  should  incorporate  the  essential  part  of  an  efficient  multigrid  solver  for 
its  corresponding  operator:  sometimes  this  essential  part  is  just  the  relaxation  part  of  that  solver,  sometimes 
this  may  even  be  an  entire  separate  multigrid  solver  applied  over  subdomains. 

4.1.  Distribution  Matrices.  For  incompressible  Navier-Stokes  equations,  an  appropriate  distribution 
matrix  corresponding  to  the  operator  L  of  (3.2)  is 

1  0  0  -a* 

0  1  0  -Sj, 

0  0  1 

0  0  0  Q, 

yielding  the  lower  triangular  operator 

Q.  0  0  0 

0  Q.  0  0 

0  0  Q,  0 

dx  dy  a^  -A  _ 

A  possible  distribution  matrix  for  the  compressible  Euler  equations  with  the  principal  linearization 
operator  L  (2.14)  is  given  by 

■  1  0  0  -^dx  0  ■ 

0  1  0  -^dy  0 

p  y 

(4.6)  M  =  0  0  1  -ia^  0 

0  0  0  Q  0 

_  0  0  0  0  1  _ 

with 


Q  0  0  0  0 

0  Q  0  0  0 

(4.7)  LM=0  0  Q  0  0. 

p(?dx  pc^dy  p(?dz  —  c^A  0 

.  py  ^dz  Q 

For  compressible  Navier-Stokes  equations,  one  of  the  factors  of  the  principal-linearization  determinant 
(2.17)  is  very  complicated.  Instead  of  devising  a  suitable  relaxation  method  for  this  scalar  factor,  one  can 
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opt  to  a  distributed  relaxation  partially  decoupling  the  linear  system  associated  with  operator  L  (2.16).  In 
particular,  the  distribution  matrix 


(4.8) 


M  = 


1  0  0  -ia*  0 

0  1  0  -ia,  0 

0  0  1  -ia^  0 

Aa,,  Xdy  Aa^  0 

p 

0  0  0  0  1 


results  in 


(4.9) 


Qr 

p 

0 

0 

0 

0 

0 

Qr 

p 

0 

0 

0 

LM  = 

0 

0 

Qr 

p 

0 

0 

Vd:, 

Vdy 

Vd, 

QQ  A+irc  ^ 

(1  —  7)kA 

_ 1 

—d 

-d, 

7 

-^A 

7P 

- 

where  V  =  pc^  +  \Q.  The  last  two  equations  remain  coupled,  requiring  a  block  2-by-2  matrix  solution  in 
relaxation.  This  distributed  relaxation  scheme  is  still  much  less  expensive  than  direct  relaxation  of  matrix 
L  requiring  solution  for  a  block  5-by-5  matrix. 


4.2.  Traditional  Factorizable  Schemes.  As  mentioned  previously,  factorizability  of  a  target  discrete 
scheme  significantly  simplifies  the  distributed  relaxation  design.  The  main  obstacle  in  this  case  to  efficient 
solution  is  that  the  discretizations  obtained  for  the  scalar  factors  in  the  discrete  determinant  are  not  always 
convenient. 


V 


V 


u  e,  p  u  e,  p  u 


V 


V 


u  e,  p  u  e,  p  u 


V 


V 


Fig.  4.1.  Staggered  arrangement  of  primitive  variables  for  Navier-Stokes  diseretization. 

4.2.1.  Staggered-Grid  Discretization  for  Navier-Stokes  Eqnations.  The  staggered-grid  dis¬ 
cretization  dating  back  to  the  mid  60’s  [27,  33,  34]  is  one  of  the  first  factorizable  discretizations  for  in¬ 
compressible  flow  equations.  Compressible  flow  discretizations  with  a  staggered  arrangement  of  variables 
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have  also  been  studied  [29,  6,  52].  A  usual  placement  of  primitive  variables  in  two  dimensions  is  depicted 
in  Figure  4.1.  With  this  staggering,  (a)  the  off-diagonal  first  derivatives  in  (3.2),  (2.14),  and  (2.16)  can  be 
approximated  as  short  central  differences;  (b)  the  second  derivatives  in  (2.16)  can  be  compositions  of  cor¬ 
responding  central  first  derivatives;  (c)  the  convection-diffusion  operators,  Q„,  can  be  approximated  by  any 
proper  discretizations  Qj).  For  discrete  factorizability,  it  is  important  to  have  the  same  discretization  for  each 
of  the  Qjy-operators  in  the  momentum  equations;  the  convection-diffusion  operators  in  other  equations  can 
be  different.  The  convection  terms  in  the  momentum  and  energy  equations  are  usually  upwind  or  upwind- 
biased;  for  simplicity,  below  we  assume  that  all  these  terms  are  the  same.  With  such  differencing,  the  discrete 
schemes  mimic  the  factorizability  property  of  the  differential  equations,  and  the  discrete  system  determi¬ 
nants  can  be  factored  as  det  L**  =  ((3[))^A^  (incompressible  Navier-Stokes)  or  det  =  {Q^Y  {Q^Q^  — 
(compressible  Euler),  where  in  three  dimensions  is  the  seven-point  /i-Laplacian,  is  an  upwind  or  up¬ 
wind  biased  discretization  of  the  convection  operators  in  the  momentum  and  energy  equations,  is  a 
convection-operator  discretization  for  the  pressure  term  in  the  fourth  equation  of  (2.14),  hence  —  c^A^ 

is  a  discrete  approximation  to  the  full-potential  operator.  The  discrete  determinant  computed  for  the  com¬ 
pressible  Navier-Stokes  equations  is  similar  to  the  differential  determinant  (2.17). 

The  discrete  distribution  matrices  follow  directly  from  the  continuous  matrices  (4.4),  (4.6),  and  (4.8). 
The  short  central  differences  are  used  for  the  approximation  of  all  the  off-diagonal  first  derivatives;  the 
convection  parts  in  the  Q-operators  are  the  same  as  those  in  the  momentum  equations.  The  resulting 
products  are  similar  to  those  for  the  continuous  problems  with  the  main  diagonals  composed  of  the 

factors  of  the  discrete  determinants. 

Distributed-relaxation  solvers  have  been  successfully  applied  to  the  staggered-grid  discretization  schemes 
for  subsonic  compressible  [52]  and  incompressible  [15,  53]  flow  problems. 

In  computing  the  Euler  system  of  equations,  the  main  disadvantages  of  the  staggered-grid  scheme  relate 
to  the  discrete  stencil  of  the  full-potential  operator.  Eor  subsonic  flow  problems,  the  downwind  differencing 
applied  for  the  term  results  in  a  full-potential-operator  stencil  that  is  somewhat  wide  (because  of  the 
term)  and  poorly  aligns  with  the  physical  (cross-stream)  anisotropies  in  approaching  the  transonic 
regime.  Eor  supersonic  flow,  where  the  problem  is  purely  hyperbolic,  the  stencil  is  not  fully  upwind  (even  if 
the  term  is  upwind  differenced)  implying  more  involved  marching  schemes. 

Recently,  a  new  approach  to  building  discretization  schemes  that  allows  any  desired  differencing  for  the 
full-potential  factor  of  the  system  determinant  without  compromising  the  scheme  factorizability  has  been 
discovered.  This  approach  is  discussed  subsequently  in  application  to  central  collocated-grid  discretizations 
(see  also  [26]),  but  it  applies  to  staggered  grids  as  well. 


4.2.2.  Collocated-Grid  Discretizations  for  the  Euler  Equations.  Another  example  of  a  factor- 
izable  scheme  is  a  collocated-grid  scheme  with  the  second-order  central  differencing  for  the  off-diagonal  first 
derivatives  in  (3.2),  (2.14),  and  (2.16).  The  convection  operators  in  the  momentum  and  energy  equations  are 
again  upwind  or  upwind  biased;  the  differencing  of  the  convection  term  applied  to  the  pressure  may  alternate 
from  downwind  (downwind-biased)  in  subsonic  mode  to  upwind  (upwind-biased)  in  supersonic  mode. 

A  typical  difficulty  associated  with  this  type  of  schemes  is  a  poor  measure  of  /i-ellipticity  in  the  discrete 
approximation  for  the  full-potential  factor  of  the  system  determinant.  To  be  more  specific,  let  us  define  the 
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collocated-grid  discretization  L**  of  the  matrix  operator  (2.14)  as 


(4.10) 


Q'^  0  0  0  ■ 

0  Q'*  0  0 

0  0  Q'*  0 

pc^dy  pc'^d^  0 

0  _ 


where  the  discrete  derivatives,  a(f,  a^ ,  a^ ,  in  all  off-diagonal  terms  are  the  wide  (with  mesh  spacing  2h)  second- 
order-accurate  central-differencing  approximations.  All  the  diagonal  terms,  except  in  the  fourth 
equation,  are  discretized  with  the  same  second-order-accurate  upwind  (or  upwind-biased)  discretization 
scheme.  In  the  subsonic  regime  (|up  =  <  c^),  the  Q^-term  is  discretized  with  a  second-order- 

accurate  downwind  (or  downwind-biased)  discretization.  The  determinant  of  the  matrix  operator  L**  is  given 

by 

(4.11)  {Q'^y  [Q'^Q'^ 

where  is  a  wide  discretization  of  the  Laplace  operator.  The  full-potential-operator  approximation 
appearing  in  the  brackets  has  two  major  drawbacks:  (1)  For  subsonic  velocities  (|u|  <  c),  the  discrete 
operator  is  not  /i-elliptic,  and  efficiency  of  any  local  relaxation  severely  degrades.  (2)  For  near-sonic  regimes 
(the  Mach  number  M  =  |u|/c  ss  1),  the  discrete  operator  stencil  does  not  reflect  the  physical  anisotropies  of 
the  differential  full-potential  operator;  the  discrete  operator  exhibits  a  very  strong  coupling  in  the  streamwise 
direction,  while  the  differential  operator  has  strong  coupling  only  in  the  cross-stream  directions. 

Several  approaches  to  cure  the  lack  of  /i-ellipticity  (mainly  in  applications  to  incompressible-flow  equa¬ 
tions)  have  been  proposed  in  the  literature  (e.g.,  [1,  13]).  Some  of  the  approaches  are  associated  with  the 
introduction  of  additional  terms  increasing  the  measure  of  /i-ellipticity  in  the  system  of  equations,  others 
propose  averaging  (filtering)  spurious  oscillations.  The  problem  of  wrong  anisotropies  in  the  full-potential- 
operator  has  not  been  sufficiently  investigated.  In  two  dimensions,  it  is  possible  to  construct  a  discretization 
that  satisfies  the  following  properties:  (1)  At  low  Mach  numbers,  the  discretization  is  dominated  by  the 
standard  (with  mesh  spacing  h)  /i-elliptic  Laplacian.  (2)  For  the  transonic  Mach  numbers,  the  discretization 
tends  to  the  optimal  discretization  [10,  11,  21]  for  the  sonic-flow  full-potential  operator.  (3)  For  supersonic 
Mach  numbers,  the  discretization  becomes  upwind  (upwind-biased)  and  can  be  solved  by  marching.  The 
problem  of  constructing  a  good  high-order  discretization  for  the  transonic  full-potential  operator  in  three 
dimensions  is  still  open. 


4.3.  Non-Factorizable  schemes.  The  majority  of  discrete  schemes  in  current  use,  especially  for  com¬ 
pressible  flow  but  also  more  recently  for  incompressible  flow,  are  based  upon  a  flux-splitting  approach.  The 
basis  of  this  approach  is  the  solution  of  the  Riemann  problem  (i.e.,  the  time  evolution  of  two  regions  of 
flow  initially  separated  by  a  diaphragm)  applied  on  a  dimension  by  dimension  basis.  This  methodology  has 
enabled  the  robust  treatment  of  flows  with  strong  shocks  and  complex  geometries.  However,  the  derived 
schemes  are  not  discretely  factorizable,  except  in  one  dimension. 

These  discrete  equations  have  always  been  solved  using  collective  relaxation  (or  pseudo-time-stepping) 
in  multidimensional  multigrid  algorithms.  A  better  efficiency  should  be  realizable  with  a  relaxation  scheme 
that  efficiently  reduces  both  the  high-frequency  and  characteristic  error  components.  Such  a  scheme  should 
combine  two  different  relaxation  methods:  (1)  A  relaxation  scheme  treating  directly  the  conservation  equa¬ 
tions  and  reducing  the  high-frequency  error  components.  (2)  A  defect-correction  (or  predictor-corrector) 
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method  with  a  factorizable  driver  (predictor)  reducing  characteristic  error  components.  This  approach  has 
not  been  tried  as  yet. 

4.4.  Boundary  Treatment.  Boundaries  and  discontinuities  introduce  some  additional  complexity  in 
distributed  relaxation.  The  determinant  of  L  M  is  usually  higher  order  than  the  determinant  of  L.  Thus,  as 
a  set  of  new  variables,  (5w  would  generally  need  additional  boundary  conditions.  In  relaxation,  because  the 
ghost  variables  can  be  added  in  the  external  part  of  the  domain,  it  is  usually  possible  to  determine  suitable 
boundary  conditions  for  Sw  that  satisfy  the  original  boundary  conditions  for  the  primitive  variables. 

Distributed  relaxation  is  applied  throughout  the  entire  computational  domain  having  the  full  effect  away 
from  boundaries  (considering  discontinuities  as  a  special  boundary  case)  in  the  regular  (smoothly  varying) 
flow  field.  The  discrete  equations  near  the  boundaries  are  usually  different  from  the  interior  equations;  the 
relaxation  equations  are  coupled  near  the  boundaries,  not  decoupled  as  they  are  in  the  interior  of  the  domain. 
Thus,  some  local  procedures  should  supplement  the  distributed  relaxation  pass.  The  coupled  near-boundary 
equations  can  be  separated  from  other  equations  and  solved  (relaxed)  with  an  appropriate  method,  such 
as  direct  solution  or  block-Newton-Kacmarcz  relaxation.  The  additional  work  will  not  seriously  affect  the 
overall  complexity  because  the  number  of  boundary  and/or  discontinuity  points  is  usually  very  small  in 
comparison  with  the  number  of  interior  points. 

Solution  (or  extensive  relaxation)  of  the  coupled  near-boundary  equations  serves  for  two  purposes.  The 
first  is  to  provide  convenient  and  reliable  boundary  conditions  for  the  distributed-relaxation  equations  in 
the  interior.  The  second  purpose  arises  because  in  the  outer  multigrid  cycle,  efficient  fine-to-coarse  grid 
restriction  of  residuals  near  the  boundaries  is  difficult  to  design;  it  depends  on  many  factors  such  as  the 
shape  of  the  boundary,  the  type  of  the  boundary  conditions,  etc.,  and  differs  from  the  residual  restriction  in 
the  interior.  A  general  way  to  avoid  efficiency  degradation  is  to  reduce  residuals  near  the  boundaries  before 
restriction  to  a  level  that  is  significantly  below  the  residual  level  characterizing  the  interior  field.  Having 
small  residuals  near  the  boundaries  makes  the  precise  form  of  the  restriction  operator  less  important. 

5.  Relaxation  of  Scalar  Factors.  Efficiency  of  the  multigrid  solvers  exploiting  factorizability  of  the 
Navier-Stokes  system  of  equations  is  determined  by  the  efficiency  of  the  relaxation  (solution)  schemes  applied 
for  scalar  factors  appearing  in  the  system  determinant.  For  uniformly  elliptic  operators  such  as  a  Laplacian, 
a  diffusion-dominated  convection-diffusion  operator,  and  a  subsonic  full-potential  operator  many  efficient  re¬ 
laxation  techniques  are  available  (see  textbooks  [6,  17,  55,  56]).  For  such  operators,  an  important  relaxation 
requirement  is  efficient  reduction  of  high-frequency  errors.  All  the  smooth  components  are  well  approximated 
on  coarse  grids  built  by  standard  (full)  coarsening;  therefore,  the  coarse-grid  correction  is  efficient  in  reduc¬ 
tion  of  smooth  errors.  For  nonelliptic  and  weakly  elliptic  factors,  e.g.,  convection,  convection-dominated 
convection-diffusion,  transonic  and  supersonic  full-potential  operators,  (smooth)  characteristic  components 
cannot  be  approximated  with  standard  multigrid  methods  [4,  10,  11,  15,  22]. 

Several  approaches  aimed  at  curing  the  characteristic-component  problem  have  been  studied  in  the 
literature.  These  approaches  fall  into  two  categories:  (1)  development  of  a  suitable  relaxation  scheme  (single¬ 
grid  method)  to  eliminate  not  only  high-frequency  error  components  but  the  characteristic  error  components 
as  well;  (2)  devising  an  adjusted  coarse-grid  operator  to  approximate  well  the  fine-grid  characteristic  error 
components. 

5.1.  Single-Grid  Methods. 

5.1.1.  Downstream  marching.  For  hyperbolic  problems,  the  simplest  first-category  method  is  down¬ 
stream  marching.  If  the  corresponding  discretization  is  a  stable  upwind  discretization  and  the  characteristic 
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field  does  not  recirculate,  then  downstream  marching  is  a  very  efficient  solver  that  yields  an  accurate  so¬ 
lution  to  a  nonlinear  hyperbolic  equation  in  just  a  few  sweeps  (a  single  downstream  sweep  provides  the 
exact  solution  to  a  linearized  problem).  The  downstream  marching  technique  was  successfully  applied  in 
solving  many  CFD  problems  associated  with  non-recirculating  flows  (see,  e.g.,  [15,  37,  45,  52,  53]).  However, 
if  a  discretization  operator  is  not  fully  upwind  (e.g.,  is  only  upwind  biased),  straightforward  downstream 
marching  is  unstable.  For  the  schemes  that  cannot  be  directly  marched,  there  are  two  possible  alternatives 
(also  of  marching  type):  defect-correction  and  predictor-corrector  methods. 

5.1.2.  Defect  Correction.  Let  us  consider  a  defect  correction  method  for  a  discretized  hyperbolic 
equation 

(^•1)  ^  f ii  ,i2  ^ 

with  specified  inflow  boundary  conditions  uo,i2  ■ 

Let  ftp, *2  be  the  current  solution  approximation.  Then  the  improved  approximation  is  calculated 
by  defect-correction  scheme  in  the  following  two  steps: 

1.  The  correction  is  calculated  by  solving  operator  with  a  right-hand  side  represented  by  the 

residual  of  (5.1)  computed  for  the  current  approximation  The  inflow  boundary  conditions  for 

V  are  initialized  with  the  zero  values. 

(^•^)  ^d^iiA2  fii,i2  T  Uii,i2‘ 

2.  The  current  approximation  is  corrected  as 

(^•H)  ^ii,i2  i  \  ^12  d2  ‘ 

The  operator  is  called  the  driver  operator.  It  is  chosen  to  be  easily  solvable  and  usually  less  accurate 
than  the  target  operator  the  latter  can  be  very  general.  If  the  iteration  converges,  steps  (5.2)  and  (5.3) 
can  be  repeated  until  the  desired  accuracy  is  reached.  Usually  the  efficiency  of  defect-correction  methods 
is  quite  satisfactory  [3,  19,  31,  52,  53],  even  though  in  principle  the  convergence  rate  of  a  defect-correction 
method  for  nonelliptic  operator  is  normally  mesh-size  dependent  [7,  23],  as  explained  below. 

In  several  papers  (e.g.,  [19,  51]),  authors  studying  the  defect-correction  method  for  nonelliptic  problems 
observed  a  slow  convergence  or  even  a  divergence  in  some  common  error  norms  for  the  initial  iterations 
and  good  asymptotic  convergence  rates  afterward.  This  behavior  is  different  from  that  observed  in  solving 
elliptic  problems  by  the  defect-correction  method,  where  the  asymptotic  convergence  rate  is  the  slowest  one. 
This  nonelliptic  feature  is  explained  by  some  properties  associated  with  the  cross-characteristic  interaction 
(e.g.,  dissipation  and/or  dispersion)  in  the  operators  involved  in  the  defect-correction  iterations.  Specifically, 
this  cross-characteristic  interaction  defines  the  penetration  distanee  (also  termed  “survival  distance”  [15]) 
of  a  characteristic  component.  The  penetration  distance  is  the  distance  from  the  inflow  boundary  within 
which  the  discrete  solution  of  the  homogeneous  problem  reasonably  approximates  the  continuous  one  (i.e., 
the  discretization  error  is  substantially  smaller  than  the  solution). 

The  penetration  distance  of  a  characteristic  component  is  roughly  proportional  to  ,  where  q  is 

the  highest  order  of  differentiation  in  the  hyperbolic  operator  under  considerations,  p  is  the  discrete-operator 
approximation  order,  to  is  the  cross-characteristic  frequency  of  the  characteristic  component,  and  h  is  the 
mesh  size.  The  ratio  of  penetration  distances  of  the  operators  and  is  an  important  factor  for  determin¬ 
ing  the  number  of  defect-correction  sweeps  required  to  reduce  the  algebraic  error  to  the  discretization-error 
level  or  to  reach  the  asymptotic  convergence  regime. 
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When  the  operators  and  have  the  same  approximation  order  {p  =  r),  efficiency  of  the  defect- 
correction  method  is  optimal  and  mesh-size  /i-independent  If  however  the  operators  and  have  different 
approximation  orders  {p  and  r,  respectively,  p  >  r),  then  efficiency  of  the  defect-correction  method  is  h- 
dependent;  i.e.,  the  maximal  number  of  sweeps  which  might  be  required  to  reduce  the  algebraic  error  to 
the  discretization-error  level  (or  to  reach  the  asymptotic  convergence  rates)  is  larger  on  fine  grids  than  on 
coarse  grids.  This  is  because  one  has  to  iterate  as  many  times  as  needed  to  attain  accuracy  up  to  the 
penetration  distance.  The  worst  (largest)  ratio  of  penetration  distances  is  obtained  for  characteristic 
components  for  which  the  penetration  distance  in  the  target  p-order  accurate  discretization  approaches  the 
depth  R  of  the  computational  domain  in  the  characteristic  direction.  It  follows  that  the  required  number  of 


iterations  is  proportional  to 


P  —  f' 
p+g 


5.1.3.  Predictor-Corrector.  One  potentially  efficient  but  yet  unexploited  method  to  overcome  grid- 
dependent  convergence  experienced  in  defect-correction  iterations  is  the  predictor-corrector  technique.  A 
detailed  look  into  the  defect-correction  iteration  reveals  that  the  computational  work  distribution  is  un¬ 
balanced:  (1)  Driver  operator  iterations  at  locations  beyond  the  penetration  distance  do  not  improve  the 
solution  approximation.  (2)  In  successive  iterations,  the  solution  approximation  near  the  inflow  boundary 
becomes  much  more  accurate  than  in  the  interior;  the  computational  efforts  spent  in  this  regions  could  be 
more  profitably  invested  at  the  accuracy  frontier. 

The  predictor-corrector  method  has  been  extensively  used  for  ordinary  and  time-dependent  differential 
equations  [18,  28];  however,  applications  for  steady-state  nonelliptic  problems  have  been  very  limited.  In 
predictor-corrector  schemes,  the  final  update  of  the  solution  at  a  particular  point  is  computed  from  the  local 
solution  of  the  target  operator.  The  solution  values  at  downstream  points  included  in  the  target-operator 
stencil  are  predicted  from  the  solution  of  the  driver  (predictor)  operator.  In  order  to  define  a  family  of 
predictor-corrector  schemes,  one  can  divide  the  computational  domain  into  several  time-like  layers;  the  first 
layer  contains  all  the  grid  points  adjacent  to  the  inflow  boundary.  Each  next  layer  is  composed  of  the  grid 
points  that  contribute  to  the  stencils  of  target  operators  defined  at  the  points  of  the  previous  layer  and  do 
not  belong  to  any  of  the  previous  layers. 

Now,  a  family  of  predictor-corrector  schemes  for  solving  the  correction  equation 


(5.4) 
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where  is  the  target  operator,  is  the  current  solution  approximation  with  residual  and 

is  the  desired  correction  function,  can  be  defined  as 
PCo'.  The  solution  of  (5.4)  is  approximated  by  solving 


(5.5) 


r/i.,  _  TDti 


This  scheme  is  identical  with  the  defect-correction  scheme. 

PCk  ■  Recursive  definition  of  the  derived  predictor-corrector  schemes  (recursion  with  respect  to  k)  can  be 
done  as  follows:  Assume  the  {j  —  l)-th  layer  have  already  been  finally  updated  in  the  current  sweep. 
Then,  to  calculate  new  values  at  the  next  j-th  layer  one  has  to  perform  the  following  three  steps: 

1.  To  predict  values  at  the  j-th  layer  with  PCu-i  scheme. 

2.  To  predict  values  at  the  {j  +  l)-th  layer  with  PCu-i  scheme. 

3.  To  update  values  at  the  j-th  layer  by  directly  relaxing  (5.4). 

The  schemes  for  A:  =  0, 1,2  have  been  tested  for  a  linearized  supersonic  full-potential  operator  [21]. 
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5.2.  Multigrid  methods. 

5.2.1.  Recirculation.  Downstream  marching  methods  are  not  viable  for  problems  with  closed  charac¬ 
teristics.  Alternative  discretization  and  solution  techniques  should  be  considered.  The  discretization  issue 
becomes  especially  important  for  flows  with  streamlines  that  do  not  start  and  end  at  boundaries,  but  con¬ 
stitute  closed  curves.  In  such  cases,  even  a  very  small  viscosity  plays  an  important  role  in  determining  the 
flow  throughout  the  domain.  The  solution  in  the  limit  of  vanishing  viscosity  depends  very  strongly  on  how 
the  viscosity  coefficients  tend  to  zero.  The  propagation  of  information  from  the  boundary  into  the  domain  is 
governed  by  the  viscous  terms  no  matter  how  small  they  may  be.  It  has  been  shown  [14]  for  both  the  scalar 
convection-diffusion  problem  and  the  incompressible  Navier-Stokes  equations  that  varying  cross-stream  nu¬ 
merical  viscosity  (caused  usually  by  varying  angles  between  the  stream  and  the  grid  lines;  e.g.,  in  standard 
upwind  and  upwind  biased  schemes)  may  prevent  convergence  to  a  physically  realizable  solution.  In  the 
most  general  case,  it  can  be  shown  that  even  isotropic  viscosity  is  not  sufficient  for  convergence  to  a  physical 
solution,  and  one  must  actually  specify  a  uniform  viscosity.  However,  for  the  homogeneous  problems  there 
are  several  indications  [14,  59]  (though  no  proof)  that  isotropy  suffices. 

To  obtain  a  discretization  scheme  that  exhibits  the  appropriate  physical-like  behavior  for  vanishing  vis¬ 
cosity,  one  must  either  add  sufficient  explicit  isotropic  viscosity  that  will  dominate  the  anisotropic  numerical 
viscosity  of  the  convection  operator,  or  else  derive  a  discrete  convection  operator  with  numerical  viscosity 
satisfying  the  condition  of  isotropy.  An  upwind  isotropic-viscosity  discretization  has  been  derived  [59] . 

One  general  approach  to  the  algebraic  solution  of  nonehiptic  equations  with  closed  characteristics  is 
to  apply  a  multigrid  method  with  an  overweighted  upwind-biased  residual  restriction.  Efficient  multigrid 
solvers  for  recirculating  convection  equation  have  already  been  demonstrated  [16,  59].  This  approach  is 
well  combined  with  the  distributed  relaxation  method  for  the  Navier-Stokes  equations,  because  within  a 
distributed  relaxation  sweep  a  multigrid  solver  with  optimal  overweighting  can  be  applied  to  a  separate 
scalar  nonehiptic  equation  with  closed  characteristics. 

Another  solution  approach  is  to  apply  some  general  techniques  to  approximate  indirectly  smooth  char¬ 
acteristic  components.  Among  helpful  techniques  are  recombination  of  iterants,  cycles  with  high  indexes, 
and  implicit  alternative-direction  relaxation.  Recombination  of  iterants  (solution  approximations  on  different 
stages  of  a  multigrid  algorithm)  at  each  grid  level  eliminates  several  (number  of  iterants  minus  one)  error  com¬ 
ponents,  those,  more  specifically,  that  are  most  prominent  in  the  residual  function.  Making  increasingly  many 
coarse-grid  iterations  per  each  fine-grid  iteration,  cycles  with  high  indexes  solve  the  characteristic-component 
problem  on  coarser  grids.  Implicit  alternative-direction  relaxation  simulates  downstream  marching  in  the  re¬ 
gions  with  open  characteristics  and  efficiently  transfers  information  in  the  regions  with  characteristics  closely 
aligned  with  the  grid.  Theoretically,  each  of  these  methods  cannot  completely  resolve  the  problem  of  poor 
coarse-grid  correction  to  the  fine-grid  smooth  characteristic  error  components.  The  problem  already  mani¬ 
fests  itself  in  two-level  algorithms  with  any  type  of  local  relaxation.  On  fine  grids  the  number  of  problematic 
error  components  may  increase,  and  many  cycles  may  be  needed  to  collect  the  necessary  number  of  the  fine- 
grid  iterants  to  exclude  all  the  troubling  error  components.  However,  it  has  been  shown  experimentally  [32] 
that  a  combination  of  implicit  alternative-direction  defect-correction  type  relaxation,  recombination  of  iter¬ 
ants  on  ah  the  levels,  and  W-cycles  can  result  in  a  relatively  efficient  multigrid  solver  for  recirculating  flow 
problems  on  practical  grids. 

5.2.2.  Full-Potential  Operator.  The  full-potential  operator  is  a  variable  type  operator,  and  its  so¬ 
lution  requires  different  procedures  in  subsonic,  transonic,  and  supersonic  regions.  In  deep  subsonic  regions, 
the  full-potential  operator  is  uniformly  elliptic  and  therefore  standard  multigrid  methods  yield  optimal 
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efficiency.  When  the  Mach  number  approaches  unity,  the  operator  becomes  increasingly  anisotropic  and,  be¬ 
cause  smooth  characteristic  error  components  cannot  be  approximated  adequately  on  coarse  grids,  classical 
multigrid  methods  severely  degrade.  In  the  deep  supersonic  regions,  the  full-potential  operator  is  uniformly 
hyperbolic  with  the  stream  direction  serving  as  the  time-like  direction.  In  this  region,  an  efficient  solver  can 
be  obtained  with  a  downstream  marching  method.  However,  downstream  marching  becomes  problematic 
for  the  Mach  number  dropping  towards  unity,  because  the  Courant  number  associated  with  this  method 
becomes  large.  Thus,  a  special  procedure  is  required  to  provide  an  efficient  solution  for  transonic  regions. 
A  possible  local  procedure  [10,  11,  20,  21]  is  based  on  piecewise  semicoarsening  and  some  rules  for  adding 
dissipation  at  the  coarse-grid  levels. 

A  similar  technique  can  be  used  to  construct  an  efficient  marching-free  multigrid  solver  for  convection- 
dominated  equations.  This  method  [22]  employs  a  colored  relaxation  scheme  and  is  very  attractive  for 
massive  parallel  computing.  A  highly  parallel  multigrid  solver  for  the  supersonic  full-potential  operator  may 
be  obtained  by  methods  similar  to  the  wave/ray  multigrid  [12]. 

6.  Analysis.  As  mentioned  previously,  it  is  important  in  attaining  optimal  efficiency  to  understand  all 
the  difficulties  that  present  themselves  in  application.  Analysis  methods  are  quite  helpful  in  this  regard,  and 
the  main  tools  are  discussed  below.  In  iterative  methods  solving  elliptic  problems,  the  main  mechanism  of 
convergence  is  damping  of  error  components.  In  solution  of  hyperbolic  scalar  equations,  there  is  another 
very  important  convergence  mechanism:  the  downstream  evolution  of  the  error  components.  In  the  presence 
of  this  additional  mechanism,  the  accuracy  first  achieved  near  the  inflow  boundary  and  is  then  propagated 
into  the  interior  of  the  domain. 

The  recognition  of  this  additional  convergence  mechanism  urges  modifications  in  the  standard  analysis 
developed  for  elliptic  problems.  Basically,  one  can  distinguish  four  types  of  analysis  applied  to  nonelliptic 
problems:  (1)  standard  linear-algebraic  matrix  analysis,  (2)  modified  zero-mode-exclusion  full-space  Fourier 
mode  analysis,  (3)  half-space  analysis  of  the  first  differential  approximations  (FDA)  [4,  57,  58],  and  (4) 
discrete  half-space  analysis.  Briefly,  the  first  differential  approximation  (also  called  modified  equation)  to  a 
difference  operator  on  a  grid  with  mesh  size  h  is  the  Taylor  expansion  of  this  difference  operator  in  terms  of 
h  truncated  to  the  first  terms  including  the  least  nonzero  power  of  h.  The  quality  of  an  analysis  applied  to 
nonelliptic  problems  is  determined  by  how  well  the  analysis  handles  the  characteristic  components. 

6.1.  Matrix  analysis.  The  most  general  and  precise  analysis  methods  are  the  linear-algebra  matrix 
analysis  methods  applied  to  the  corresponding  linearized  problem.  This  analysis  considers  the  difference 
operators  without  assumptions  about  the  solution  and  boundary  conditions.  It  can  be  applied  to  variable- 
coefficient  problems  as  well.  This  analysis  was  found  very  useful  for  analyzing  one-dimensional  problems. 
However,  the  enormous  computational  complexity  of  this  analysis  makes  it  not  viable  for  multidimensional 
problems.  Although,  the  analysis  complexity  can  be  reduced  considerably  by  assuming  Fourier  modes  in  two 
of  the  three  spatial  directions. 

6.2.  Modified  full-space  Fourier  mode  analysis.  The  modified  full-space  Fourier  mode  analysis  is 
a  modification  of  the  standard  full-space  Fourier  mode  analysis  excluding  from  the  consideration  all  the  zero 
modes  (the  modes  with  vanishing  symbols)  [56].  It  is  the  simplest  and  most  popular  type  of  analysis  (e.g., 
see  applications  in  [19,  31]).  This  analysis  estimates  only  the  amplification  (damping)  factor.  Its  inherent 
disadvantage  is  the  inability  to  take  the  influence  of  the  inflow  boundary  into  account.  This  explains  its 
failure  in  describing  the  downstream  error  evolution.  However,  the  modified  full-space  analysis  can  also  be 
useful  for  analyzing  the  effect  of  forcing  terms. 
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6.3.  FDA  half-space  analysis.  The  FDA  half-space  analysis  is  a  relatively  simple  and  efficient  tool 
for  analyzing  the  effect  of  the  inflow  boundary.  Examples  of  applications  of  this  analysis  are  available  [4,  15, 
57,  58].  The  first  differential  approximations  are  considered  on  a  half-space  including  an  inflow  boundary. 
The  boundary  conditions  are  represented  by  one  Fourier  mode  at  a  time.  The  FDA  analysis  provides 
a  good  qualitative  description  of  the  downstream  error  evolution.  This  analysis  focuses  on  characteristic 
components  and,  therefore,  considers  homogeneous  problems.  Note  that  a  combination  of  the  FDA  analysis 
with  the  modified  full-space  analysis  can  provide  a  good  insight  for  nonhomogeneous  problems  as  well.  The 
disadvantages  of  this  analysis  are  the  inability  to  provide  quantitative  estimates,  to  analyze  the  effect  of 
different  boundary  condition  discretizations,  and  to  address  the  asymptotic  convergence  rate. 

6.4.  Discrete  half-space  analysis.  The  discrete  half-space  analysis  [10,  23]  considers  the  discretiza¬ 
tions  in  their  exact  form  rather  than  their  differential  approximation,  while  the  boundary  data  are  rep¬ 
resented  by  a  Fourier  component.  This  analysis  translates  the  original  multidimensional  problem  into  a 
one-dimensional  discrete  problem,  where  the  frequencies  of  the  boundary  Fourier  components  are  considered 
as  parameters.  To  regularize  the  half-space  problem,  the  solution  is  not  allowed  to  grow  faster  than  a  poly¬ 
nomial  function.  This  tool  is  very  accurate;  it  can  be  used  to  explain  in  detail  many  phenomena  observed 
in  solving  nonelliptic  equations  and  provides  a  close  prediction  of  the  actual  solution  behavior. 

The  one-dimensional  solution  obtained  in  the  discrete  half-space  analysis  has  two  different  representation 
forms:  (1)  away  from  the  boundary,  the  solution  is  defined  as  a  linear  combination  of  a  finite  number  of 
analytical  components;  this  region  is  called  the  analytical  representation  region;  (2)  in  the  region  adjacent  to 
the  inflow  boundary,  the  solution  is  defined  pointwise;  this  zone  is  referred  to  as  the  pointwise  representation 
region.  With  each  further  iteration  described  by  the  analysis,  the  pointwise  representation  region  penetrates 
by  a  finite  number  of  mesh  sizes  into  the  interior.  By  using  these  representations,  the  computational  com¬ 
plexity  of  the  analysis  becomes  much  less  than  that  associated  with  the  one-dimensional  matrix  analysis.  In 
the  asymptotic  regime,  when  the  pointwise  representation  zone  covers  all  the  domain,  this  analysis  becomes 
a  discrete  one-dimensional  matrix  analysis  of  the  multidimensional  problem. 

The  discrete  half-space  analysis  provides  a  quantitative  description  of  the  approximate  solution;  it  pre¬ 
dicts  the  convergence  rate  for  each  iteration  and  the  asymptotic  convergence  rate.  It  can  be  easily  adjusted 
to  analyze  the  global  effect  of  any  local  discretization  of  the  inflow  boundary  conditions.  This  adjustment 
can  be  done  just  by  widening  the  initial  pointwise  representation  region  at  the  inflow  boundary.  If  neces¬ 
sary,  the  analysis  can  take  into  account  the  influence  of  the  discretized  outflow  boundary  conditions  as  well. 
Generally,  this  discrete  half-space  analysis  treats  completely  both  mechanisms  of  convergence,  damping  and 
downstream  evolution  of  errors,  associated  with  nonelliptic  problem  solvers. 

7.  Summary  of  Recent  Progress. 

7.1.  Pressure-Equation  Discretization.  The  original  pressure-equation  formulation  [45]  has  been 
extended  to  general  coordinates  and  implemented  for  lifting  airfoils  in  inviscid  flow  [36,  37,  39]  and  viscous 
flow  [47].  The  results  for  viscous  flow  over  a  lifting  airfoil  at  low  Reynolds  number  are  shown  in  Figure  7.1. 
An  alternating  line-implicit  Gauss-Seidel  relaxation  is  used  to  treat  the  mesh  anisotropy  that  generally 
occurs  in  resolving  viscous  boundary  layers  on  stretched  grids.  The  computed  pressure  distributions  are 
nearly  indistinguishable  from  each  other  on  the  finer  grids.  The  convergence  rate  actually  becomes  better 
as  grids  are  refined  and  more  levels  in  the  FAS  cycle  are  used;  the  16x8  grid  is  always  the  coarsest  grid  in 
the  multigrid  computations.  The  convergence  rates  are  comparable  to  the  rates  obtained  for  fully  elliptic 
problems. 
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Karman-Trefftz  Airfoil,  Re^  =  200,  a  =  2® 
Surface  Pressure 


(a)  Pressure  distribution  on  the  upper  and  lower  surfaces  for  a  sequence  of  grids. 


Karman-Trefftz  Airfoil,  Re^  =  200,  a  =  2° 


(b)  Residual  convergence  with  an  FMG  method  using  10  FAS  cycles  on  each  grid  from  the  coarsest  16x8  to 
the  finest  256x128  (squaresia;— momentum;  trianglesiy— momentum;  circles ipressure  equation). 

Fig.  7.1.  Computational  results  for  the  ineompressible  viseous  flow  over  a  lifting  Karman-Trefftz  airfoil  at  Re  =  200  and 
a  =  2°. 

7.2.  Staggered-Grid  Factorizable  Discretization.  The  first  TME  solver  applying  the  distributed 
relaxation  approach  for  solution  to  an  entering  flow  problem  for  the  incompressible  Navier-Stokes  equations 
was  developed  using  a  staggered-grid  formulation  [15].  This  formulation  was  extended  to  the  compressible 
NS  equations  and  fast  convergence  rates  were  demonstrated  [52]  for  several  viscous  model  problems.  This 
latter  work  was  the  first  experience  with  distributed  relaxation  in  the  computation  of  compressible  viscous 
flows  in  which  a  2  x  2  block  was  relaxed  simultaneously  in  line-implicit  Gauss-Seidel  relaxation.  The  coupling 
of  boundary  and  interior  relaxation  was  not  optimally  treated  at  the  time.  A  more  complete  study  on  TME 
for  the  incompressible  equations  at  high-Reynolds-number  conditions  has  been  recently  performed  [53].  In 
all  calculations,  a  staggered  arrangement  of  variables  on  Cartesian  grids  has  been  used.  With  distributed 
relaxation,  the  system  of  equations  has  been  decomposed  (i.e.,  factored)  everywhere,  except  near  boundaries 
where  the  equations  remained  coupled.  The  results  of  the  calculation  are  shown  in  Eigure  7.2  for  the  viscous 
flow  over  a  finite  fiat  plate.  The  convergence  of  residuals  and  the  algebraic-to-discretization  errors  in  drag 
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Cycles 


(a)  Algebraic-to-discretization  errors  in  drag,  Cd- 


(b)  L2— norm  of  the  residual. 

Fig.  7.2.  Convergence  of  errors  in  an  FMG  cycle  using  5  FAS  cycles  on  each  grid  from  the  coarsest  13x7  to  the  finest 
193x97  for  the  incompressible  viscous  flow  over  a  flat  plate  at  Re  =  10,  000. 

are  shown  versus  multigrid  cycles.  The  residual  convergence  rate  is  about  the  same  as  for  the  underlying 
Laplacian  factor.  The  FMG  solver  with  just  one  FAS  multigrid  cycle  per  grid  level  and  a  total  computational 
work  equivalent  to  about  10  target-grid  residual  evaluations  converged  the  drag  to  the  discretization  accuracy. 

7.3.  New  Factorizable  Collocated-Grid  Discretizations.  Recently,  a  new  multidimensional  fac- 
torizable  scheme  for  the  Euler  equations  has  been  developed  [42]  for  Cartesian  coordinates  and  extended 
through  generalized  coordinates  to  external  lifting  flows  around  airfoils  with  both  subcritical  and  supercriti¬ 
cal  freestream  Mach  numbers  [38,  35].  The  starting  point  for  the  scheme  is  the  first-order  discretization  of  the 
flux-difference  splitting  scheme  [40] .  Correction  terms  are  added  in  the  form  of  mixed  derivatives  to  make  the 
scheme  both  second-order  accurate  and  discretely  factorizable.  The  resulting  scheme  is  second-order  accurate 
and  compact  in  comparison  to  other  scheme;  it  is  the  first  fiux-difference-splitting  scheme  that  is  discretely 
factorizable  in  multiple  dimensions.  Discrete  factorizability  is  achieved  by  using  some  non-standard  wide 
approximations  for  spatial  derivatives  to  ensure  that  the  identities 
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(a)  Surface  pressure  (Cp)  on  a  sequence  of  three  grids. 


(b)  Convergence  of  residual  and  forces  (lift  and  drag)  in  an  FMG  cycle  with  20  FAS  cycles  on  the  finest 
257  X  257  mesh  and  10  FAS  cycles  on  each  coarser  mesh. 

Fig.  7.3.  Computational  results  for  compressible  Euler  flow  over  a  lifting  Kdrmdn-Trefftz  airfoil;  Moo  =  0.70;  a  =  1° . 


^xx^yy  —  ^xy^xy: 

^xx^y  —  ^xy^x: 

^yy^x  —  ^xy^y 

are  satisfied  on  the  discrete  level.  The  determinant  of  the  resulting  scheme  is  composed  of  an  upwind 
differenced  convection  factor  and  an  /i-elliptic  approximation  for  the  full-potential  factor.  The  distributed 
relaxation  is  possible  by  using  a  left  and  right  distribution  matrix,  although  this  has  not  been  applied  as  yet. 

In  numerical  tests  performed  for  this  scheme,  the  multigrid  solver  employed  symmetric  point  collec¬ 
tive  Gauss-Seidel  relaxation.  Computations  for  subsonic  and  transonic  channel  flows  with  essentially  grid- 
independent  convergence  rates  have  been  presented  [38] .  Grid-independent  convergence  rates  have  also  been 
attained  for  a  flow  with  stagnation  points  [35].  The  convergence  rates  observed  in  experimenting  with 
subsonic  flow  over  a  lifting  airfoil  were  quite  fast  (about  0.3  per  multigrid  V-cycle)  and  only  slightly  grid 
dependent.  The  rates  somewhat  deteriorate  in  transonic/supersonic  computations,  emphasizing  the  need 
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for  distributed  relaxation.  The  scheme  applies  at  low  Mach  numbers  although  it  has  yet  to  be  extended 
to  viscous  flows.  Multigrid  results  for  the  transonic  flow  over  a  lifting  Karman-Trefftz  airfoil  with  a  shock 
are  shown  in  Figure  7.3.  The  pressure  distribution  shows  the  weak  shock  that  is  captured  by  the  scheme. 
The  residual  convergence  indicates  some  deterioration  of  the  rate  on  the  finer  grids  but  the  lift  and  drag 
coefficients  are  converged  to  below  discretization  error  levels  in  only  a  few  cycles. 

Another  approach  to  building  factorizable  schemes  with  suitable  discretizations  for  scalar  factors  has 
been  explored  in  papers  of  the  second  and  third  authors  [24,  26,  25].  The  approach  is  based  on  a  collocated- 
grid  scheme  with  a  mechanism  that  allows  one  to  improve  the  /i-ellipticity  measure  by  obtaining  any  desired 
discretizations  for  the  full-potential  factor  of  the  system  determinant  without  compromising  the  discrete 
factorizability.  Also,  the  distribution  matrices  follow  directly  from  the  discrete  forms  for  M  presented  earlier. 
The  same  approach  can  be  applied  for  incompressible-flow  problems  and  to  staggered-grid  discretizations  as 
well. 

The  starting  point  is  the  discretization  (4.10).  The  way  proposed  to  improve  the  discrete  full-potential 
operator  is  to  change  the  discretization  of  to  +  A'^.  Then  the  discrete  full-potential  operator  is 
changed  to 

(7.1)  Q'^A'^  +Q'^Q'^ 

where  A^  =  “  {Q^Q^  —  c^A^^),  and  is  a  desired  approximation  for  the  full- 

potential  factor.  In  smooth  regions,  A^  is  second-order  small  (proportional  to  /i^),  hence  the  overall  second- 
order  discretization  accuracy  is  not  compromised.  The  operator  (<3'*)”^  is  a  nonlocal  operator  and  its 
introduction  can  be  effected  through  a  new  auxiliary  variable  and  a  new  discrete  equation  =  V^p^. 

Thus,  the  corrected  discrete  approximation  to  (2.14)  is  defined  as 

'  Q'*  0  0  0  0  ' 

0  0  0  0 

(7  2)  0  0  0  0 

0  0  0  Q'*  -V^  0 

pc^dy  pc^d^  1  0 

.  0  0  Q\ 

The  corresponding  distribution  matrix,  M**,  for  distributed  relaxation  is  defined  as 

1  0  0  0  0  ' 

0  10  0  -iaf)  0 

p  y 

0  0  10  -iaj  0 

P  Z 

0  0  0  1  v'^  o’ 

0  0  0  0  Q'^  0 

0  0  0  0  0  1 

so  that  the  resulting  matrix  L**  M**  becomes  lower  triangular  as 

g'‘  0  0  0  0  o' 

0  g'‘  0  0  0  0 

0  0  g'‘  0  0  0 

0  0  0  g'‘  0  o' 

pAdx  pc^dy  pAd'p  1  .7^'“  0 

Ad’p  Ad’p  0  g'* 

'y  J-  'y  y  'y  Z  J 
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The  scheme  as  defined  above  is  valid  for  nonconservative  flows.  A  version  to  be  used  for  distributed 
relaxation  of  conservative  equations  has  also  been  designed  [24] . 

Numerical  tests  have  only  been  performed  as  yet  for  a  quasi-one-dimensional  subsonic  flow  in  a  con¬ 
vergent/divergent  channel.  The  accuracy  was  comparable  to  other  schemes.  With  proper  treatment  of  the 
distributed-relaxation  equations  in  the  regions  adjacent  to  the  boundaries,  the  convergence  of  the  multigrid 
solver  with  a  V-cycle  and  two  relaxation  sweeps  per  level  is  identical  with  the  convergence  of  a  similar 
multigrid  solver  for  the  discrete  full-potential  operator. 

8.  Concluding  Remarks.  Fundamentals  and  recent  advances  towards  the  development  of  TME  solvers 
for  fluid  simulations  have  been  presented.  Accurate  discrete  approximations  to  the  solution  of  the  differential 
equations  are  obtained  with  FMG  methods  through  fast  reduction  of  algebraic  errors  below  the  discretization 
error  level  on  each  mesh.  Strategies  to  attain  TME  for  general  fluid  systems  by  exploiting  factorizability 
of  the  governing  differential  equations  are  reviewed.  These  strategies  include  a  reformulation  of  the  target 
differential  equations  and  a  distributed  relaxation  approach  applied  to  the  original  equations.  New  discretiza¬ 
tions  and  computations  demonstrating  this  methodology  for  inviscid  and  viscous  flow  simulations  have  been 
presented. 


REFERENCES 

[1]  S.  W.  Armfield,  Ellipticity,  Accuracy,  and  Convergence  of  the  Discrete  Navier-Stokes  Equations,  J. 

Comput.  Phys.,  114  (1994),  pp.  176-184. 

[2]  D.  Bai  and  a.  Brandt,  Local  mesh  refinement  multilevel  techniques,  SIAM  J.  Sci.  Stat.  Comput.,  8 

(1987),  pp.  109-135. 

[3]  K.  Bohmer,  P.  W.  Hemker,  and  H.  J.  Stetter,  The  defect  correction  approach,  in  Defect  Cor¬ 

rection  Methods,  K.  Bohmer  and  H.  J.  Stetter,  eds..  Comp.  Suppl.  5,  Wien,  New  York,  1984, 
Springer-Verlag,  pp.  1-32. 

[4]  A.  Brandt,  Multigrid  solvers  for  non-elliptic  and  singular-perturbation  steady-state  problems,  (unpub¬ 

lished).  The  Weizmann  Institute  of  Science,  Rehovot,  Israel,  December  1981. 

[5]  - ,  Guide  to  multigrid  development,  in  Multigrid  Methods,  W.  Hackbusch  and  U.  Trottenberg,  eds.. 

Lecture  Notes  in  Math.  960,  Springer-Verlag,  Berlin,  1982. 

[6]  - ,  Multigrid  techniques:  1984  guide  with  applications  to  fluid  dynamics,  in  Lecture  Notes  for  the 

Computational  Fluid  Dynamics,  Lecture  Series  at  the  Von-Karman  Institute  for  Fluid  Dynamics, 
The  Weizmann  Institute  of  Science,  Rehovot,  Israel,  1984.  ISBN-3-88457-081-1,  GMD-Studien  Nr. 
85,  Available  from  GMD-AIW,  Postfach  1316,  D-53731,  St.  Augustin  1,  Germany.  Also  available 
from  Secretary,  Department  of  Mathematics,  University  of  Colorado  at  Denver,  CO  80204-5300. 

[7]  - ,  The  Weizmann  Institute  of  Science  research  in  multilevel  computations:  1988  Report,  inProc.  4th 

Copper  Mountain  Conf.  on  Multigrid  Methods,  J.  Mandel  and  et  ah,  eds.,  SIAM,  1989,  pp.  13-53. 

[8]  - ,  Barriers  to  achieving  textbook  multigrid  efficiency  in  CED,  ICASE  Interim  Report  32,  NASA 

CR-1998-207647,  April  1998. 

[9]  - ,  Appendix  C:  Recent  developments  in  multigrid  efficiency  in  computational  fluid  dynamics,  in 

Multigrid,  Academic  Press,  London,  2000,  pp.  573-589.  Ulrich  Trottenberg,  C.  W.  Oosterlee,  and 
A.  Schuler. 

[10]  A.  Brandt  and  B.  Diskin,  Multigrid  solvers  for  the  non-aligned  sonic  flow:  The  constant  coefficient 


24 


case,  Computers  and  Fluids,  28  (1999),  pp.  511-549.  Also  Gauss  Center  Report  WI/GC-8  at  The 
Weizmann  Institute  of  Science,  Rehovot,  Israel,  October  1997. 

[11]  - ,  Multigrid  solvers  for  nonaligned  sonic  flows,  SIAM  J.  Sci.  Comp.,  21  (2000),  pp.  473-501. 

[12]  A.  Brandt  and  I.  Livshits,  Wave-ray  multigrid  methods  for  standing  wave  equations.  Electronic 

Trans.  Num.  An.,  6  (1997),  pp.  162-181. 

[13]  A.  Brandt  and  S.  Ta’asan,  Multigrid  solutions  to  quasi- elliptic  schemes,  in  Progress  and  Supercom¬ 

puting  in  Computational  Fluid  Dynamics,  E.  M.  Murman  and  S.  S.  Abarbanel,  eds.,  Boston,  MA, 
1985,  Birkhauser,  pp.  143-154. 

[14]  A.  Brandt  and  I.  Yavneh,  Inadequacy  of  first-order  upwind  difference  schemes  for  some  recirculating 

flows,  J.  Comput.  Phys.,  93  (1991),  pp.  128-143. 

[15]  - ,  On  multigrid  solution  of  high-Reynolds  incompressible  entering  flow,  J.  Comput.  Phys.,  101 

(1992),  pp.  151-164. 

[16]  - ,  Accelerated  multigrid  convergence  and  high-Reynolds  recirculating  flows,  SIAM  J.  Sci.  Comp.,  14 

(1993),  pp.  607-626. 

[17]  W.  L.  Briggs,  S.  F.  McCormick,  and  V.  E. Henson,  Multigrid  Tutorial,  2nd  edition,  SIAM,  USA, 

2000. 

[18]  K.  Burrage,  Parallel  and  Sequential  Methods  for  Ordinary  Differential  Equations,  Claredon  Press, 

Oxford,  1995. 

[19]  J.  A.  Desideri  and  P.  W.  Hemker,  Convergence  analysis  of  the  defect- correction  iteration  for  hy¬ 

perbolic  problems,  SIAM  J.  Sci.  Comp.,  16  (1995),  pp.  88-118. 

[20]  B.  Diskin,  Multigrid  algorithm  with  conditional  coarsening  for  the  non-aligned  sonic  flow,  Eiectronic 

Trans.  Num.  An.,  6  (1997),  pp.  106-119.  Aiso  in  Proceedings  of  the  Eighth  Copper  Mountain  Con¬ 
ference  on  Multigrid  Methods,  1997. 

[21]  - ,  Efficient  Multigrid  Solvers  for  the  Linearaized  Transonic  Eull  Potential  Equation,  PhD  thesis. 

The  Weizmann  Institute  of  Science,  1998. 

[22]  - ,  Efficient  multigrid  methods  for  solving  upwind-biased  discretizations  of  the  convection  equation, 

Appiied  Mathematics  and  Computation,  123  (2001),  pp.  343-379.  (Aiso  ICASE  Report  99-25,  NASA 
CR-1999-209355). 

[23]  B.  Diskin  and  J.  L.  Thomas,  Half-space  analysis  of  the  defect- correction  method  for  Eromm  dis¬ 

cretization  of  convection,  SIAM  J.  Sci.  Comp.,  22  (2000),  pp.  633-655. 

[24]  B.  Diskin  and  J.  L.  Thomas,  Distributed  relaxation  for  conservative  discretizations,  AIAA  Paper 

2001-2571,  15th  AIAA  CFD  Conference,  Anaheim,  CA,  June  2001. 

[25]  B.  Diskin  and  J.  L.  Thomas,  Analysis  of  boundary  conditions  for  factorizable  discretizations  of  the 

Euler  equations,  ICASE  Report  2002-13,  NASA  CR-2002-211648,  May  2002. 

[26]  - ,  New  factorizable  discretizations  for  the  Euler  equations,  ICASE  Report  2002-6,  NASA  CR-2002- 

211456,  April  2002.  Submitted  to  SIAM  J.  Sci.  Comp. 

[27]  F.  H.  Harlow  and  J.  E.  Welch,  Numerical  calculations  of  time- dependent  viscous  incompressible 

flow  of  fluid  with  free  surface.  Physics  of  Fiuids,  8  (1965),  pp.  2182-2189. 

[28]  C.  Hirsch,  Numerical  computation  of  internal  and  external  flows.  Vol.l,  Eundumentals  of  numerical 

discretization,  A  Wiiey-Interscience  pubiication,  John  Wiiey  &  Sons,  Inc.,  605  Third  Avenue,  New 
York,  NY  10158-0012,  USA,  1988. 

[29]  K.  C.  Karki  and  S.  V.  Patankar,  Pressure  based  calculation  procedure  for  viscous  flows  at  all  speeds 

in  arbitrary  configurations,  AIAA  Journal,  27  (1989),  pp.  1167-1174. 


25 


[30]  R.  S.  Montero,  I.  M.  Llorente,  and  M.  D.  Salas,  Robust  muHigrid  algorithm  for  the  Navier- 

Stokes  equations,  J.  Comput.  Phys.,  173  (2001),  pp.  412-432. 

[31]  C.  W.  OoSTERLEE,  F.  J.  Gaspar,  T.  Washio,  AND  R.  WiENANDS,  Multigrid  line  smoothers  for 

higher  order  upwind  diseretizations  of  eonveetion-dominated  problems,  J.  Comput.  Phys.,  139  (1998), 
pp.  274-307. 

[32]  C.  W.  OoSTERLEE  AND  T.  Washio,  Krylov  subspaee  aeeeleration  of  nonlinear  multigrid  with  appliea- 

tion  to  reeireulating  flows,  SIAM  J.  Scient.  Comp.,  21(5)  (2000),  pp.  1670-1690. 

[33]  S.  V.  Patankar,  Numerieal  Heat  Transfer  and  Fluid  Flow,  Hemisphere  Publishing  Co. /McGraw-Hill 

Co.,  New  York,  1980. 

[34]  R.  Peyret  and  T.  D.  Taylor,  Computational  Methods  for  Fluid  Flow,  Springer  Verlag,  New  York, 

1983. 

[35]  T.  W.  Roberts,  The  development  of  a  faetorizable  multigrid  algorithm  for  subsonie  and  transonie  flow, 

AIAA  Paper  2001-2572,  15th  AIAA  CFD  meeting,  Anaheim,  CA,  June  2001. 

[36]  T.  W.  Roberts,  D.  Sidilkover,  and  R.  C.  Swanson,  Textbook  multigrid  effieieney  for  the  steady 

Euler  equations,  AIAA  Paper  97-1949,  13th  AIAA  CFD  meeting,  Snowmass  Village,  CO,  June-July 
1997. 

[37]  - ,  An  algorithm  for  ideal  multigrid  eonvergenee  for  the  steady  Euler  equations.  Computers  and 

Fluids,  28  (1999),  pp.  427-442. 

[38]  T.  W.  Roberts,  D.  Sidilkover,  and  J.  L.  Thomas,  Multigrid  relaxation  of  a  faetorizable  eonser- 

vative  diseretization  of  the  eompressible  Euler  equations,  June  2000.  AIAA  Paper  2000-2252. 

[39]  T.  W.  Roberts,  D.  Sidilkover,  and  S.  V.  Tsynkov,  On  the  eombined  performanee  of  nonloeal 

artifieial  boundary  eonditions  with  the  new  generation  of  advaneed  multigrid  flow  solvers.  Computers 
and  Fluids,  31  (2002),  pp.  269-308. 

[40]  P.  L.  Roe,  Charaeteristie-based  sehemes  for  the  Euler  equations,  in  Ann.  Rev.  Fluid  Mech.,  vol.  18, 

1986,  pp.  337-365. 

[41]  R.  Rubinstein,  C.  L.  Rumsey,  M.  D.  Salas,  and  J.  L.  Thomas,  Turbulenee  modeling  workshop, 

ICASE  Interim  Report  37,  NASA  CR-2001-210841,  March  2001. 

[42]  D.  Sidilkover,  Faetorizable  seheme  for  the  equation  of  fluid  flow,  ICASE  Report  99-20,  NASA  CR- 

1999-209345,  June  1999. 

[43]  - ,  Some  approaehes  toward  eonstrueting  optimally  effieient  multigrid  solvers  for  the  inviseid  flow 

equations.  Computers  and  Fluids,  28  (1999),  pp.  551-571. 

[44]  D.  Sidilkover,  Faetorizable  upwind  sehemes:  The  triangular  unstruetured  grid  formulation,  AIAA 

Paper  2001-2575,  15th  AIAA  CFD  meeting,  Anaheim,  CA,  June  2001. 

[45]  D.  Sidilkover  and  U.  Asher,  A  multigrid  solver  for  the  steady-state  Navier-Stokes  equations  using 

the  pres  sure- Poisson  formulation,  Matematica  Aplicada  e  Computational,  14  (1995),  pp.  21-35. 

[46]  K.  Stuben  and  U.  Trottenberg,  Multigrid  methods:  Fundamental  algorithms,  model  problem  anal¬ 

ysis  and  applieation,  in  Multigrid  Methods,  W.  Hackbusch  and  U.  Trottenberg,  eds..  Lecture  Notes 
in  Math.  960,  Springer- Verlag,  Berlin,  1982,  pp.  1-176. 

[47]  R.  C.  Swanson,  Towards  optimal  multigrid  effieieney  for  the  Navier-Stokes  equations,  AIAA  Paper 

2001-2574,  15th  AIAA  CFD  Conference,  Anaheim,  CA,  June  2001. 

[48]  S.  Ta’asan,  Canonieal  forms  of  multidimensional  steady  inviseid  flows,  ICASE  Report  93-34,  NASA 

CR-191488,  1993. 

[49]  - ,  CanonieaTvariables  multigrid  method  for  steady-state  Euler  equations,  ICASE  Report  94-14, 


26 


NASA  CR-194888,  1994. 

[50]  - ,  Essentially  optimal  multigrid  method  for  steady  state  Euler  equations,  AIAA  Paper  95-0209,  33rd 

Aerospace  Sciences  Meeting  and  Exhibit,  January  1995. 

[51]  J.  L.  Thomas,  D.  L.  Bonhaus,  W.  K.  Anderson,  C.  L.  Rumsey,  and  R.  T.  Biedron,  An  0(nrn^) 

plane  solver  for  the  eompressihle  Navier-Stokes  equations,  AIAA  Paper  99-0785,  37th  Aerospace 
Sciences  Meeting  &  Exhibit,  Reno,  NV,  January  1999. 

[52]  J.  L.  Thomas,  B.  Diskin,  and  A.  Brandt,  Distributed  relaxation  multigrid  and  defeet  eorreetion 

applied  to  the  eompressihle  Navier-Stokes  equations,  AIAA  Paper  99-3334,  14th  Computational  Eluid 
Dynamics  Conference,  Norfolk,  VA,  July  1999. 

[53]  - ,  Textbook  multigrid  effieieney  for  the  ineompressible  Navier-Stokes  equations:  High  Reynolds  num¬ 

ber  wakes  and  boundary  layers.  Computers  and  Eluids,  30  (2001),  pp.  853-874.  Also  ICASE  Report 
99-51  (NASA  CR-1999-209831),  December  1999. 

[54]  J.  L.  Thomas,  B.  Diskin,  A.  Brandt,  and  J.  C.  South,  General  framework  for  aehieving  textbook 

multigrid  effieieney:  Quasi-l-d  Euler  example,  in  Erontiers  of  Computational  Eluid  Dynamics  — 
2002,  D.  A.  Caughey  and  M.  M.  Hafez,  eds..  World  Scientific  Publishing  Company,  Singapore,  2002, 
pp.  61-79.  Also  ICASE  Report  2000-30,  NASA/CR-2000-210320. 

[55]  U.  Trottenberg,  C.  W.  Oosterlee,  and  A.  Schuler,  Multigrid,  Academic  Press,  London,  2000. 

[56]  P.  Wesseling,  An  introduetion  to  multigrid  methods.  Pure  and  Applied  Mathematics,  John  Wiley  & 

Sons,  Chichester,  1992. 

[57]  N.  Yanenko  and  Y.I.Shokin,  Eirst  differential  approximation  method  and  approximate  viseosity  of 

differenee  sehemes.  High-speed  eomputing  in  fluid  dynamies.  The  Physics  of  Eluids,  Supplement  II, 
New  York,  1969. 

[58]  - ,  On  the  first  differential  approximations  of  differenee  sehemes  for  hyperbolie  systems  of  equations, 

Siberian  Mathematical  Journal,  10  (1969),  pp.  1174-1188. 

[59]  I.  Yavneh,  C.  Venner,  and  A.  Brandt,  East  multigrid  solution  of  the  adveetion  problem  with  elosed 

eharaeteristies,  SIAM  J.  Sci.  Comp.,  19  (1998),  pp.  111-125. 


27 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
0MB  No.  0704-0188 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources, 
gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this 
collection  of  information,  including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson 
Davis  Highway,  Suite  1204,  Arlington,  VA  22202-4302,  and  to  the  Office  of  Management  and  Budget,  Paperwork  Reduction  Project  (0704-0188),  Washington,  DC  20503. 

1.  AGENCY  USE  ONLYfieai/e  blank) 

2.  REPORT  DATE 

May  2002 

3.  REPORT  TYPE  AND  DATES  COVERED 

Contractor  Report 

4.  TITLE  AND  SUBTITLE 

Recent  advances  in  achieving  textbook  multigrid  efficiency  for 
computational  fluid  dynamics  simulations 

5.  FUNDING  NUMBERS 

C  NASl-97046 

WU  505-90-52-01 

6.  AUTHOR(S) 

Achi  Brandt,  Boris  Diskin,  and  James  L.  Thomas 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

ICASE 

Mail  Stop  132C 

NASA  Langley  Research  Center 

Hampton,  VA  23681-2199 

8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

ICASE  Report  No.  2002-16 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

National  Aeronautics  and  Space  Administration 

Langley  Research  Center 

Hampton,  VA  23681-2199 

10.  SPONSORING/MONITORING 

AGENCY  REPORT  NUMBER 

NASA/CR-2002-211656 

ICASE  Report  No.  2002-16 

11.  SUPPLEMENTARY  NOTES 

Langley  Technical  Monitor:  Dennis  M.  Bushnell 

Final  Report 

AIAA  Conference  and  Annual  Review  of  Fluid  Mechanics 

12a.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Unclassified-Unlimited 

Subject  Category  64 

Distribution:  Nonstandard 

Availability:  NASA-CASI  (301)  621-0390 

12b.  DISTRIBUTION  CODE 

13.  ABSTRACT  (Maximum  200  words) 

Recent  advances  in  achieving  textbook  multigrid  efficiency  for  fluid  simulations  are  presented.  Textbook  multigrid 
efficiency  is  defined  as  attaining  the  solution  to  the  governing  system  of  equations  in  a  computational  work  which 
is  a  small  multiple  of  the  operation  counts  associated  with  discretizing  the  system.  Strategies  are  reviewed  to 
attain  this  efficiency  by  exploiting  the  factorizability  properties  inherent  to  a  range  of  fluid  simulations,  including 
the  compressible  Navier-Stokes  equations.  Factorizability  is  used  to  separate  the  elliptic  and  hyperbolic  factors 
contributing  to  the  target  system;  each  of  the  factors  can  then  be  treated  individually  and  optimally.  Boundary 
regions  and  discontinuities  are  addressed  with  separate  (local)  treatments.  New  formulations  and  recent  calculations 
demonstrating  the  attainment  of  textbook  efficiency  for  aerodynamic  simulations  are  shown. 

14.  SUBJECT  TERMS 

textboook  multigrid  efficiency,  factorizable  systems  of  differential  equations, 
principal  linearization,  factorizable  discretizations,  distributed  relaxation 

15.  NUMBER  OF  PAGES 

32 

16.  PRICE  CODE 

A03 

17.  SECURITY  CLASSIFICATION 

OF  REPORT 

Unclassified 

18.  SECURITY  CLASSIFICATION 
OF  THIS  PAGE 

Unclassified 

19.  SECURITY  CLASSIFICATION 
OF  ABSTRACT 

20.  LIMITATION 

OF  ABSTRACT 

NSN  7540-01-280-5500  Standard  Form  298(Rev.  2-89) 

Prescribed  by  ANSI  Std.  Z39-18 


298-102 


